TrackingTools/test/MuonErrorMatrixAnalyzer.cc

0001 #include "RecoMuon/TrackingTools/test/MuonErrorMatrixAnalyzer.h"
0002
0003 #include "DataFormats/GeometrySurface/interface/Plane.h"
0004 #include "TrackingTools/PatternTools/interface/TSCPBuilderNoMaterial.h"
0005 //#include <TrackingTools/GeomPropagators/interface/AnalyticalImpactPointExtrapolator.h>
0006 #include "TrackingTools/TrajectoryParametrization/interface/CurvilinearTrajectoryParameters.h"
0007
0008 //#include <TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h>
0009 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0010
0011 #include "DataFormats/TrackReco/interface/Track.h"
0012 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0013 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0014
0015 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0016
0017 #include <cmath>
0018 #include "DataFormats/Math/interface/deltaPhi.h"
0019
0020 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
0021
0022 #include "MagneticField/Engine/interface/MagneticField.h"
0023 #include "SimDataFormats/Associations/interface/TrackToTrackingParticleAssociator.h"
0024
0025 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
0026 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
0027
0028 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0029 #include "FWCore/ServiceRegistry/interface/Service.h"
0030
0031 #include "TH1F.h"
0032 #include "TH2F.h"
0033 #include "TF1.h"
0034 #include "TROOT.h"
0035 #include "TList.h"
0036 #include "TKey.h"
0037
0038 //
0039 // constructors and destructor
0040 //
0041 MuonErrorMatrixAnalyzer::MuonErrorMatrixAnalyzer(const edm::ParameterSet& iConfig) {
0042   theCategory = "MuonErrorMatrixAnalyzer";
0043   //now do what ever initialization is needed
0044   theTrackLabel = iConfig.getParameter<edm::InputTag>("trackLabel");
0045
0046   trackingParticleLabel = iConfig.getParameter<edm::InputTag>("trackingParticleLabel");
0047
0048   //adding this for associator
0049   theAssocLabel = iConfig.getParameter<std::string>("associatorName");
0050
0051   theErrorMatrixStore_Reported_pset = iConfig.getParameter<edm::ParameterSet>("errorMatrix_Reported_pset");
0052
0053   theErrorMatrixStore_Residual_pset = iConfig.getParameter<edm::ParameterSet>("errorMatrix_Residual_pset");
0054   theErrorMatrixStore_Pull_pset = iConfig.getParameter<edm::ParameterSet>("errorMatrix_Pull_pset");
0055
0056   thePlotFileName = iConfig.getParameter<std::string>("plotFileName");
0057
0058   theRadius = iConfig.getParameter<double>("radius");
0059   if (theRadius != 0) {
0060     GlobalPoint O(0, 0, 0);
0061     Surface::RotationType R;
0062     refRSurface = Cylinder::build(theRadius, O, R);
0063     thePropagatorName = iConfig.getParameter<std::string>("propagatorName");
0064     thePropagatorToken = esConsumes(edm::ESInputTag("", thePropagatorName));
0065     theZ = iConfig.getParameter<double>("z");
0066     if (theZ != 0) {
0067       //plane can only be specified if R is specified
0068       GlobalPoint Opoz(0, 0, theZ);
0069       GlobalPoint Oneg(0, 0, -theZ);
0070       refZSurface[1] = Plane::build(Opoz, R);
0071       refZSurface[0] = Plane::build(Oneg, R);
0072     }
0073   }
0074
0075   theGaussianPullFitRange = iConfig.getUntrackedParameter<double>("gaussianPullFitRange", 2.0);
0076 }
0077
0078 MuonErrorMatrixAnalyzer::~MuonErrorMatrixAnalyzer() {
0079   // do anything here that needs to be done at desctruction time
0080   // (e.g. close files, deallocate resources etc.)
0081 }
0082
0083 //
0084 // member functions
0085 //
0086
0087 // ------------ method called to for each event  ------------
0088 void MuonErrorMatrixAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0089   if (theErrorMatrixStore_Reported) {
0090     analyze_from_errormatrix(iEvent, iSetup);
0091   }
0092
0093   if (theErrorMatrixStore_Residual || theErrorMatrixStore_Pull) {
0094     analyze_from_pull(iEvent, iSetup);
0095   }
0096 }
0097
0098 FreeTrajectoryState MuonErrorMatrixAnalyzer::refLocusState(const FreeTrajectoryState& fts) {
0099   if (theRadius == 0) {
0100     GlobalPoint vtx(0, 0, 0);
0101     TSCPBuilderNoMaterial tscpBuilder;
0102     FreeTrajectoryState PCAstate = tscpBuilder(fts, vtx).theState();
0103     return PCAstate;
0104   } else {
0105     //go to the cylinder surface, along momentum
0106     TrajectoryStateOnSurface onRef = thePropagator->propagate(fts, *refRSurface);
0107
0108     if (!onRef.isValid()) {
0109       edm::LogError(theCategory) << " cannot propagate to cylinder of radius: " << theRadius;
0110       //try out the plane if specified
0111       if (theZ != 0) {
0112         onRef = thePropagator->propagate(fts, *(refZSurface[(fts.momentum().z() > 0)]));
0113         if (!onRef.isValid()) {
0114           edm::LogError(theCategory) << " cannot propagate to the plane of Z: "
0115                                      << ((fts.momentum().z() > 0) ? "+" : "-") << theZ << " either.";
0116           return FreeTrajectoryState();
0117         }  //invalid state
0118       }    //z plane is set
0119       else {
0120         return FreeTrajectoryState();
0121       }
0122     }  //invalid state
0123     else if (fabs(onRef.globalPosition().z()) > theZ && theZ != 0) {
0124       //try out the plane
0125       onRef = thePropagator->propagate(fts, *(refZSurface[(fts.momentum().z() > 0)]));
0126       if (!onRef.isValid()) {
0127         edm::LogError(theCategory) << " cannot propagate to the plane of Z: " << ((fts.momentum().z() > 0) ? "+" : "-")
0128                                    << theZ << " even though cylinder z indicates it should.";
0129       }  //invalid state
0130     }    //z further than the planes
0131
0132     LogDebug(theCategory) << "reference state is:\n" << onRef;
0133
0134     return (*onRef.freeState());
0135   }  // R=0
0136 }
0137
0138 void MuonErrorMatrixAnalyzer::analyze_from_errormatrix(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0139   using namespace edm;
0140
0141   if (theRadius != 0) {
0142     //get a propagator
0143     thePropagator = iSetup.getHandle(thePropagatorToken);
0144   }
0145
0146   //get the mag field
0147   theField = iSetup.getHandle(theFieldToken);
0148
0149   //open a collection of track
0150   edm::Handle<reco::TrackCollection> tracks;
0151   iEvent.getByLabel(theTrackLabel, tracks);
0152
0153   //loop over them
0154   for (unsigned int it = 0; it != tracks->size(); ++it) {
0155     //   take their initial free state
0156     FreeTrajectoryState PCAstate = trajectoryStateTransform::initialFreeState((*tracks)[it], theField.product());
0157     if (PCAstate.position().mag() == 0) {
0158       edm::LogError(theCategory) << "invalid state from track initial state. skipping.\n" << PCAstate;
0159       continue;
0160     }
0161
0162     FreeTrajectoryState trackRefState = refLocusState(PCAstate);
0163     if (trackRefState.position().mag() == 0)
0164       continue;
0165
0166     AlgebraicSymMatrix55 errorMatrix = trackRefState.curvilinearError().matrix();
0167
0168     double pT = trackRefState.momentum().perp();
0169     double eta = fabs(trackRefState.momentum().eta());
0170     double phi = trackRefState.momentum().phi();
0171
0172     LogDebug(theCategory) << "error matrix:\n" << errorMatrix << "\n state: \n" << trackRefState;
0173
0174     //fill the TProfile3D
0175     for (int i = 0; i != 5; ++i) {
0176       for (int j = i; j != 5; ++j) {
0177         //get the profile plot to fill
0178         TProfile3D* ij = theErrorMatrixStore_Reported->get(i, j);
0179         if (!ij) {
0180           edm::LogError(theCategory) << "cannot get profile " << i << " " << j;
0181           continue;
0182         }
0183
0184         //get sigma squared or correlation factor
0185         double value = MuonErrorMatrix::Term(errorMatrix, i, j);
0186         ij->Fill(pT, eta, phi, value);
0187       }
0188     }
0189   }
0190 }
0191
0192 void MuonErrorMatrixAnalyzer::analyze_from_pull(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0193   using namespace edm;
0194
0195   //get the mag field
0196   theField = iSetup.getHandle(theFieldToken);
0197
0198   //open a collection of track
0199   edm::Handle<View<reco::Track> > tracks;
0200   iEvent.getByLabel(theTrackLabel, tracks);
0201
0202   //open a collection of Tracking particle
0203   edm::Handle<TrackingParticleCollection> TPtracks;
0204   iEvent.getByLabel(trackingParticleLabel, TPtracks);
0205
0206   //get the associator
0207   edm::Handle<reco::TrackToTrackingParticleAssociator> associator;
0208   iEvent.getByLabel(theAssocLabel, associator);
0209
0210   //associate
0211   reco::RecoToSimCollection recSimColl = associator->associateRecoToSim(tracks, TPtracks);
0212
0213   LogDebug(theCategory) << "I have found: " << recSimColl.size() << " associations in total.";
0214
0215   int it = 0;
0216   for (reco::RecoToSimCollection::const_iterator RtSit = recSimColl.begin(); RtSit != recSimColl.end(); ++RtSit) {
0217     //what am I loop over ?
0218     //the line below this has a problem
0219     //    const reco::TrackRef & track = RtSit->key;
0220     const std::vector<std::pair<TrackingParticleRef, double> >& tp = RtSit->val;
0221
0222     //what do I want to get from those
0223     FreeTrajectoryState sim_fts;
0224     FreeTrajectoryState trackPCAstate;
0225
0226     LogDebug(theCategory) << "I have found: " << tp.size() << " tracking particle associated with this reco::Track.";
0227
0228     if (tp.size() != 0) {
0229       //take the match with best quality
0230       std::vector<std::pair<TrackingParticleRef, double> >::const_iterator vector_iterator = tp.begin();
0231       const std::pair<TrackingParticleRef, double>& pair_in_vector = *vector_iterator;
0232       //what am I looking at ?
0233       const TrackingParticleRef& trp = pair_in_vector.first;
0234       //the line below this has a syntax error
0235       //double quality & = pair_in_vector.second;
0236       //     const TrackingParticle & best_associated_trackingparticle = trp.product();
0237
0238       //     work on the TrackingParticle
0239       //get reference point and momentum to make fts
0240       GlobalPoint position_sim_fts(trp->vx(), trp->vy(), trp->vz());
0241       GlobalVector momentum_sim_fts(trp->px(), trp->py(), trp->pz());
0242       int charge_sim = (int)trp->charge();
0243       GlobalTrajectoryParameters par_sim(position_sim_fts, momentum_sim_fts, charge_sim, theField.product());
0244
0245       sim_fts = FreeTrajectoryState(par_sim);
0246
0247       //       work on the reco::Track
0248       trackPCAstate = trajectoryStateTransform::initialFreeState((*tracks)[it], theField.product());
0249       if (trackPCAstate.position().mag() == 0) {
0250         edm::LogError(theCategory) << "invalid state from track initial state. skipping.\n" << trackPCAstate;
0251         continue;
0252       }
0253
0254       //get then both at the reference locus
0255       FreeTrajectoryState simRefState = refLocusState(sim_fts);
0256       FreeTrajectoryState trackRefState = refLocusState(trackPCAstate);
0257
0258       if (simRefState.position().mag() == 0 || trackRefState.position().mag() == 0)
0259         continue;
0260
0261       GlobalVector momentum_sim = simRefState.momentum();
0262       GlobalPoint point_sim = simRefState.position();
0263
0264       GlobalVector momentum_track = trackRefState.momentum();
0265       GlobalPoint point_track = trackRefState.position();
0266
0267       //conversion from global to curvinlinear parameters for both reco and sim track
0268       CurvilinearTrajectoryParameters sim_CTP(point_sim, momentum_sim, simRefState.charge());
0269       CurvilinearTrajectoryParameters track_CTP(point_track, momentum_track, trackRefState.charge());
0270
0271       //These are the parameters for the CTP point
0272       AlgebraicVector5 sim_CTP_vector = sim_CTP.vector();
0273       AlgebraicVector5 track_CTP_vector = track_CTP.vector();
0274       const AlgebraicSymMatrix55& track_error = trackRefState.curvilinearError().matrix();
0275
0276       double pT_sim = simRefState.momentum().perp();
0277
0278       //get the momentum, eta, phi here
0279       double pT = trackRefState.momentum().perp();
0280       double eta = fabs(trackRefState.momentum().eta());
0281       double phi = trackRefState.momentum().phi();
0282
0283       LogDebug(theCategory) << "The sim pT for this association is: " << pT_sim << " GeV"
0284                             << "sim state: \n"
0285                             << simRefState << "The track pT for this association is: " << pT << " GeV"
0286                             << "reco state: \n"
0287                             << trackRefState;
0288
0289       //once have the momentum,eta,phi choose what bin it should go in
0290       //how do i get the correct bin
0291       double diff_i = 0, diff_j = 0;
0292       for (unsigned int i = 0; i != 5; ++i) {
0293         //do these if statements because if the parameter is phi it is not simply reco minus sim
0294         if (i != 2) {
0295           diff_i = sim_CTP_vector[i] - track_CTP_vector[i];
0296         } else {
0297           diff_i = deltaPhi(sim_CTP_vector[i], track_CTP_vector[i]);
0298         }
0299
0300         for (unsigned int j = i; j < 5; ++j) {
0301           //do these if statements because if the parameter is phi it is not simply reco minus sim
0302           if (j != 2) {
0303             diff_j = sim_CTP_vector[j] - track_CTP_vector[j];
0304           } else {
0305             diff_j = deltaPhi(sim_CTP_vector[j], track_CTP_vector[j]);
0306           }
0307
0308           //filling residual histograms
0309           if (theErrorMatrixStore_Residual) {
0310             unsigned int iH = theErrorMatrixStore_Residual->index(i, j);
0311             TProfile3D* ij = theErrorMatrixStore_Residual->get(i, j);
0312             if (!ij) {
0313               edm::LogError(theCategory) << i << " " << j << " not vali indexes. TProfile3D not found.";
0314               continue;
0315             }
0316
0317             int iPt = theErrorMatrixStore_Residual->findBin(ij->GetXaxis(), pT) - 1;    // between 0 and maxbin-1;
0318             int iEta = theErrorMatrixStore_Residual->findBin(ij->GetYaxis(), eta) - 1;  // between 0 and maxbin-1;
0319             int iPhi = theErrorMatrixStore_Residual->findBin(ij->GetZaxis(), phi) - 1;  // between 0 and maxbin-1;
0320
0321             TH1ptr& theH = (theHist_array_residual[iH])[index(ij, iPt, iEta, iPhi)];
0322
0323             if (i == j) {
0324               LogDebug(theCategory) << "filling for: " << i;
0325               ((TH1*)theH)->Fill(diff_i);
0326             } else {
0327               LogDebug(theCategory) << "filling for: " << i << " " << j;
0328               ((TH2*)theH)->Fill(diff_i, diff_j);
0329             }
0330           }  //filling residual histograms
0331
0332           //filling pull histograms
0333           if (theErrorMatrixStore_Pull) {
0334             unsigned int iH = theErrorMatrixStore_Pull->index(i, j);
0335             TProfile3D* ij = theErrorMatrixStore_Pull->get(i, j);
0336             if (!ij) {
0337               edm::LogError(theCategory) << i << " " << j << " not vali indexes. TProfile3D not found.";
0338               continue;
0339             }
0340
0341             int iPt = theErrorMatrixStore_Pull->findBin(ij->GetXaxis(), pT) - 1;    // between 0 and maxbin-1;
0342             int iEta = theErrorMatrixStore_Pull->findBin(ij->GetYaxis(), eta) - 1;  // between 0 and maxbin-1;
0343             int iPhi = theErrorMatrixStore_Pull->findBin(ij->GetZaxis(), phi) - 1;  // between 0 and maxbin-1;
0344
0345             TH1ptr& theH = (theHist_array_pull[iH])[index(ij, iPt, iEta, iPhi)];
0346
0347             if (i == j) {
0348               LogDebug(theCategory) << "filling pulls for: " << i;
0349               ((TH1*)theH)->Fill(diff_i / sqrt(track_error(i, i)));
0350             } else {
0351               LogDebug(theCategory) << "filling pulls (2D) for: " << i << " " << j;
0352               ((TH2*)theH)->Fill(diff_i / sqrt(track_error(i, i)), diff_j / sqrt(track_error(j, j)));
0353             }
0354           }  // filling the pull histograms
0355         }
0356       }  //loop over all the 15 terms of the 5x5 matrix
0357
0358     }  //end of if (tp.size()!=0)
0359     it++;
0360   }  //end of loop over the map
0361 }
0362
0363 // ------------ method called once each job just before starting event loop  ------------
0364
0365 void MuonErrorMatrixAnalyzer::beginJob() {
0366   if (theErrorMatrixStore_Reported_pset.empty()) {
0367     theErrorMatrixStore_Reported = 0;
0368   } else {
0369     //create the error matrix provider, saying that you want to construct the things from here
0370     theErrorMatrixStore_Reported = new MuonErrorMatrix(theErrorMatrixStore_Reported_pset);
0371   }
0372
0373   if (theErrorMatrixStore_Residual_pset.empty()) {
0374     theErrorMatrixStore_Residual = 0;
0375   } else {
0376     //create the error matrix provider for the alternative method
0377     theErrorMatrixStore_Residual = new MuonErrorMatrix(theErrorMatrixStore_Residual_pset);
0378   }
0379
0380   if (theErrorMatrixStore_Pull_pset.empty()) {
0381     theErrorMatrixStore_Pull = 0;
0382   } else {
0383     //create the error matrix provider for the alternative method
0384     theErrorMatrixStore_Pull = new MuonErrorMatrix(theErrorMatrixStore_Pull_pset);
0385   }
0386
0387   if (thePlotFileName == "") {
0388     thePlotFile = 0;
0389     //    thePlotDir=0;
0390     gROOT->cd();
0391   } else {
0392     edm::Service<TFileService> fs;
0393     //    thePlotDir = new TFileDirectory(fs->mkdir(thePlotFileName.c_str()));
0394     thePlotFile = TFile::Open(thePlotFileName.c_str(), "recreate");
0395     thePlotFile->cd();
0396     theBookKeeping = new TList;
0397   }
0398
0399   //book histograms in this section
0400   //you will get them back from their name
0401
0402   //Since the bin size is only specified in ErrorMatrix.cc must get bin size from the TProfile3D
0403   //need to choose which TProfile to get bin content info from.
0404
0405   //create the 15 histograms, 5 of them TH1F, 10 of them TH2F
0406   if (theErrorMatrixStore_Residual) {
0407     for (unsigned int i = 0; i != 5; ++i) {
0408       for (unsigned int j = i; j < 5; ++j) {
0409         unsigned int iH = theErrorMatrixStore_Residual->index(i, j);
0410         TProfile3D* ij = theErrorMatrixStore_Residual->get(i, j);
0411         if (!ij) {
0412           edm::LogError(theCategory) << i << " " << j << " not valid indexes. TProfile3D not found.";
0413           continue;
0414         }
0415         unsigned int pTBin = (ij->GetNbinsX());
0416         unsigned int etaBin = (ij->GetNbinsY());
0417         unsigned int phiBin = (ij->GetNbinsZ());
0418         //allocate memory for all the histograms, for the given matrix term
0419         theHist_array_residual[iH] = new TH1ptr[maxIndex(ij)];
0420
0421         //book the histograms now
0422         for (unsigned int iPt = 0; iPt < pTBin; ++iPt) {
0423           for (unsigned int iEta = 0; iEta < etaBin; iEta++) {
0424             for (unsigned int iPhi = 0; iPhi < phiBin; iPhi++) {
0425               TString hname = Form("%s_%i_%i_%i", ij->GetName(), iPt, iEta, iPhi);
0426               LogDebug(theCategory) << "preparing for: " << hname << "\n"
0427                                     << "trying to access at: " << index(ij, iPt, iEta, iPhi)
0428                                     << "while maxIndex is: " << maxIndex(ij);
0429               TH1ptr& theH = (theHist_array_residual[iH])[index(ij, iPt, iEta, iPhi)];
0430
0431               const int bin[5] = {100, 100, 100, 5000, 100};
0432               const double min[5] = {-0.05, -0.05, -0.1, -0.005, -10.};
0433               const double max[5] = {0.05, 0.05, 0.1, 0.005, 10.};
0434
0435               if (i == j) {
0436                 //        TString v(ErrorMatrix::vars[i]);
0437                 TString htitle(Form("%s Pt:[%.1f,%.1f] Eta:[%.1f,%.1f] Phi:[%.1f,%.1f]",
0438                                     MuonErrorMatrix::vars[i].Data(),
0439                                     ij->GetXaxis()->GetBinLowEdge(iPt + 1),
0440                                     ij->GetXaxis()->GetBinUpEdge(iPt + 1),
0441                                     ij->GetYaxis()->GetBinLowEdge(iEta + 1),
0442                                     ij->GetYaxis()->GetBinUpEdge(iEta + 1),
0443                                     ij->GetZaxis()->GetBinLowEdge(iPhi + 1),
0444                                     ij->GetZaxis()->GetBinUpEdge(iPhi + 1)));
0445                 //diagonal term
0446                 thePlotFile->cd();
0447                 theH = new TH1F(hname, htitle, bin[i], min[i], max[i]);
0448                 //        theH = thePlotDir->make<TH1F>(hname,htitle,bin[i],min[i],max[i]);
0449                 theBookKeeping->Add(theH);
0450                 theH->SetXTitle("#Delta_{reco-gen}(" + MuonErrorMatrix::vars[i] + ")");
0451
0452                 LogDebug(theCategory) << "creating a TH1F " << hname << " at: " << theH;
0453               } else {
0454                 TString htitle(Form("%s Pt:[%.1f,%.1f] Eta:[%.1f,%.1f] Phi:[%.1f,%.1f]",
0455                                     ij->GetTitle(),
0456                                     ij->GetXaxis()->GetBinLowEdge(iPt + 1),
0457                                     ij->GetXaxis()->GetBinUpEdge(iPt + 1),
0458                                     ij->GetYaxis()->GetBinLowEdge(iEta + 1),
0459                                     ij->GetYaxis()->GetBinUpEdge(iEta + 1),
0460                                     ij->GetZaxis()->GetBinLowEdge(iPhi + 1),
0461                                     ij->GetZaxis()->GetBinUpEdge(iPhi + 1)));
0462                 thePlotFile->cd();
0463                 theH = new TH2F(hname, htitle, bin[i], min[i], max[i], 100, min[j], max[j]);
0464                 //        theH = thePlotDir->make<TH2F>(hname,htitle,bin[i],min[i],max[i],100,min[j],max[j]);
0465                 theBookKeeping->Add(theH);
0466                 theH->SetXTitle("#Delta_{reco-gen}(" + MuonErrorMatrix::vars[i] + ")");
0467                 theH->SetYTitle("#Delta_{reco-gen}(" + MuonErrorMatrix::vars[j] + ")");
0468
0469                 LogDebug(theCategory) << "creating a TH2 " << hname << " at: " << theH;
0470               }
0471             }
0472           }
0473         }  //end of loop over the pt,eta,phi
0474       }
0475     }
0476   }
0477
0478   //create the 15 histograms, 5 of them TH1F, 10 of them TH2F
0479   if (theErrorMatrixStore_Pull) {
0480     for (unsigned int i = 0; i != 5; ++i) {
0481       for (unsigned int j = i; j < 5; ++j) {
0482         unsigned int iH = theErrorMatrixStore_Pull->index(i, j);
0483         TProfile3D* ij = theErrorMatrixStore_Pull->get(i, j);
0484         if (!ij) {
0485           edm::LogError(theCategory) << i << " " << j << " not valid indexes. TProfile3D not found.";
0486           continue;
0487         }
0488         unsigned int pTBin = (ij->GetNbinsX());
0489         unsigned int etaBin = (ij->GetNbinsY());
0490         unsigned int phiBin = (ij->GetNbinsZ());
0491         //allocate memory for all the histograms, for the given matrix term
0492         theHist_array_pull[iH] = new TH1ptr[maxIndex(ij)];
0493
0494         //book the histograms now
0495         for (unsigned int iPt = 0; iPt < pTBin; ++iPt) {
0496           for (unsigned int iEta = 0; iEta < etaBin; iEta++) {
0497             for (unsigned int iPhi = 0; iPhi < phiBin; iPhi++) {
0498               TString hname = Form("%s_p_%i_%i_%i", ij->GetName(), iPt, iEta, iPhi);
0499               LogDebug(theCategory) << "preparing for: " << hname << "\n"
0500                                     << "trying to access at: " << index(ij, iPt, iEta, iPhi)
0501                                     << "while maxIndex is: " << maxIndex(ij);
0502               TH1ptr& theH = (theHist_array_pull[iH])[index(ij, iPt, iEta, iPhi)];
0503
0504               const int bin[5] = {200, 200, 200, 200, 200};
0505               const double min[5] = {-10, -10, -10, -10, -10};
0506               const double max[5] = {10, 10, 10, 10, 10};
0507
0508               if (i == j) {
0509                 //              TString v(ErrorMatrix::vars[i]);
0510                 TString htitle(Form("%s Pt:[%.1f,%.1f] Eta:[%.1f,%.1f] Phi:[%.1f,%.1f]",
0511                                     MuonErrorMatrix::vars[i].Data(),
0512                                     ij->GetXaxis()->GetBinLowEdge(iPt + 1),
0513                                     ij->GetXaxis()->GetBinUpEdge(iPt + 1),
0514                                     ij->GetYaxis()->GetBinLowEdge(iEta + 1),
0515                                     ij->GetYaxis()->GetBinUpEdge(iEta + 1),
0516                                     ij->GetZaxis()->GetBinLowEdge(iPhi + 1),
0517                                     ij->GetZaxis()->GetBinUpEdge(iPhi + 1)));
0518                 //diagonal term
0519                 thePlotFile->cd();
0520                 theH = new TH1F(hname, htitle, bin[i], min[i], max[i]);
0521                 //        theH = thePlotDir->make<TH1F>(hname,htitle,bin[i],min[i],max[i]);
0522                 theBookKeeping->Add(theH);
0523                 theH->SetXTitle("#Delta_{reco-gen}/#sigma(" + MuonErrorMatrix::vars[i] + ")");
0524
0525                 LogDebug(theCategory) << "creating a TH1F " << hname << " at: " << theH;
0526               } else {
0527                 TString htitle(Form("%s Pt:[%.1f,%.1f] Eta:[%.1f,%.1f] Phi:[%.1f,%.1f]",
0528                                     ij->GetTitle(),
0529                                     ij->GetXaxis()->GetBinLowEdge(iPt + 1),
0530                                     ij->GetXaxis()->GetBinUpEdge(iPt + 1),
0531                                     ij->GetYaxis()->GetBinLowEdge(iEta + 1),
0532                                     ij->GetYaxis()->GetBinUpEdge(iEta + 1),
0533                                     ij->GetZaxis()->GetBinLowEdge(iPhi + 1),
0534                                     ij->GetZaxis()->GetBinUpEdge(iPhi + 1)));
0535                 thePlotFile->cd();
0536                 theH = new TH2F(hname, htitle, bin[i], min[i], max[i], 100, min[j], max[j]);
0537                 //        theH = thePlotDir->make<TH2F>(hname,htitle,bin[i],min[i],max[i],100,min[j],max[j]);
0538                 theBookKeeping->Add(theH);
0539                 theH->SetXTitle("#Delta_{reco-gen}/#sigma(" + MuonErrorMatrix::vars[i] + ")");
0540                 theH->SetYTitle("#Delta_{reco-gen}/#sigma(" + MuonErrorMatrix::vars[j] + ")");
0541
0542                 LogDebug(theCategory) << "creating a TH2 " << hname << " at: " << theH;
0543               }
0544             }
0545           }
0546         }  //end of loop over the pt,eta,phi
0547       }
0548     }
0549   }
0550 }
0551
0552 #include <TF2.h>
0553
0554 MuonErrorMatrixAnalyzer::extractRes MuonErrorMatrixAnalyzer::extract(TH2* h2) {
0555   extractRes res;
0556
0557   //FIXME. no fitting procedure by default
0558   if (h2->GetEntries() < 1000000) {
0559     LogDebug(theCategory) << "basic method. not enough entries (" << h2->GetEntries() << ")";
0560     res.corr = h2->GetCorrelationFactor();
0561     res.x = h2->GetRMS(1);
0562     res.y = h2->GetRMS(2);
0563     return res;
0564   }
0565
0566   //make a copy while rebinning
0567   int nX = std::max(1, h2->GetNbinsX() / 40);
0568   int nY = std::max(1, h2->GetNbinsY() / 40);
0569   LogDebug(theCategory) << "rebinning: " << h2->GetName() << " by: " << nX << " " << nY;
0570   TH2* h2r = h2->Rebin2D(nX, nY, "hnew");
0571
0572   TString fname(h2->GetName());
0573   fname += +"_fit_f2";
0574   TF2* f2 = new TF2(
0575       fname, "[0]*exp(-0.5*(((x-[1])/[2])**2+((y-[3])/[4])**2 -2*[5]*(x-[1])*(y-[3])/([4]*[2])))", -10, 10, -10, 10);
0576   f2->SetParameters(h2->Integral(), 0, h2->GetRMS(1), 0, h2->GetRMS(2), h2->GetCorrelationFactor());
0577   f2->FixParameter(1, 0);
0578   f2->FixParameter(3, 0);
0579
0580   if (fabs(h2->GetCorrelationFactor()) < 0.001) {
0581     LogDebug(theCategory) << "correlations neglected: " << h2->GetCorrelationFactor() << " for: " << h2->GetName();
0582     f2->FixParameter(5, 0);
0583   } else {
0584     f2->ReleaseParameter(5);
0585   }
0586
0587   f2->SetParLimits(2, 0, 10 * h2->GetRMS(1));
0588   f2->SetParLimits(4, 0, 10 * h2->GetRMS(2));
0589
0590   h2r->Fit(fname, "nqL");
0591
0592   res.corr = f2->GetParameter(5);
0593   res.x = f2->GetParameter(2);
0594   res.y = f2->GetParameter(4);
0595
0596   LogDebug(theCategory) << "\n variable: " << h2->GetXaxis()->GetTitle() << "\n RMS= " << h2->GetRMS(1)
0597                         << "\n fit= " << f2->GetParameter(2) << "\n variable: " << h2->GetYaxis()->GetTitle()
0598                         << "\n RMS= " << h2->GetRMS(2) << "\n fit= " << f2->GetParameter(4) << "\n correlation"
0599                         << "\n correlation factor= " << h2->GetCorrelationFactor()
0600                         << "\n fit=                " << f2->GetParameter(5);
0601
0602   f2->Delete();
0603   h2r->Delete();
0604
0605   return res;
0606 }
0607
0608 // ------------ method called once each job just after ending the event loop  ------------
0609 void MuonErrorMatrixAnalyzer::endJob() {
0610   LogDebug(theCategory) << "endJob begins.";
0611   //  std::cout<<"endJob of MuonErrorMatrixAnalyzer"<<std::endl;
0612   //evaluate the histograms to find the correlation factors and sigmas
0613
0614   if (theErrorMatrixStore_Reported) {
0615     //close the error matrix method object
0616     theErrorMatrixStore_Reported->close();
0617   }
0618
0619   //write the file with all the plots in it.
0620   TFile* thePlotFile = 0;
0621   if (thePlotFileName != "") {
0622     //    std::cout<<"trying to write in: "<<thePlotFileName<<std::endl;
0623
0624     thePlotFile = TFile::Open(thePlotFileName.c_str(), "recreate");
0625     thePlotFile->cd();
0626     TListIter iter(theBookKeeping);
0627     //    std::cout<<"number of objects to write: "<<theBookKeeping->GetSize()<<std::endl;
0628     TObject* o = 0;
0629     while ((o = iter.Next())) {
0630       //    std::cout<<"writing: "<<o->GetName()<<" in file: "<<thePlotFile->GetName()<<std::endl;
0631       o->Write();
0632     }
0633     //thePlotFile->Write();
0634   }
0635
0636   if (theErrorMatrixStore_Residual) {
0637     //extract the rms and correlation factor from the residual
0638     for (unsigned int i = 0; i != 5; ++i) {
0639       for (unsigned int j = i; j < 5; ++j) {
0640         unsigned int iH = theErrorMatrixStore_Residual->index(i, j);
0641         TProfile3D* ij = theErrorMatrixStore_Residual->get(i, j);
0642         TProfile3D* ii = theErrorMatrixStore_Residual->get(i, i);
0643         TProfile3D* jj = theErrorMatrixStore_Residual->get(j, j);
0644         if (!ij) {
0645           edm::LogError(theCategory) << i << " " << j << " not valid indexes. TProfile3D not found.";
0646           continue;
0647         }
0648         if (!ii) {
0649           edm::LogError(theCategory) << i << " " << i << " not valid indexes. TProfile3D not found.";
0650           continue;
0651         }
0652         if (!jj) {
0653           edm::LogError(theCategory) << j << " " << j << " not valid indexes. TProfile3D not found.";
0654           continue;
0655         }
0656
0657         unsigned int pTBin = (ij->GetNbinsX());
0658         unsigned int etaBin = (ij->GetNbinsY());
0659         unsigned int phiBin = (ij->GetNbinsZ());
0660
0661         //analyze the histograms now
0662         for (unsigned int iPt = 0; iPt < pTBin; ++iPt) {
0663           for (unsigned int iEta = 0; iEta < etaBin; iEta++) {
0664             for (unsigned int iPhi = 0; iPhi < phiBin; iPhi++) {
0665               double pt = ij->GetXaxis()->GetBinCenter(iPt + 1);
0666               double eta = ij->GetYaxis()->GetBinCenter(iEta + 1);
0667               double phi = ij->GetZaxis()->GetBinCenter(iPhi + 1);
0668
0669               TH1ptr& theH = (theHist_array_residual[iH])[index(ij, iPt, iEta, iPhi)];
0670               LogDebug(theCategory) << "extracting for: " << pt << " " << eta << " " << phi
0671                                     << "at index: " << index(ij, iPt, iEta, iPhi)
0672                                     << "while maxIndex is: " << maxIndex(ij) << "\n named: " << theH->GetName();
0673
0674               //FIXME. not using the i=j plots (TH1F)
0675               if (i != j) {
0676                 extractRes r = extract((TH2*)theH);
0677                 ii->Fill(pt, eta, phi, r.x);
0678                 jj->Fill(pt, eta, phi, r.y);
0679                 ij->Fill(pt, eta, phi, r.corr);
0680                 LogDebug(theCategory) << "for: " << theH->GetName() << " rho is: " << r.corr << " sigma x= " << r.x
0681                                       << " sigma y= " << r.y;
0682               }
0683
0684               //please free the memory !!!
0685               LogDebug(theCategory) << "freing memory of: " << theH->GetName();
0686               theH->Delete();
0687               theH = 0;
0688             }
0689           }
0690         }  //end of loop over the pt,eta,phi
0691       }
0692     }
0693
0694     theErrorMatrixStore_Residual->close();
0695   }
0696
0697   if (theErrorMatrixStore_Pull) {
0698     // extract the scale factors from the pull distribution
0699     TF1* f = new TF1("fit_for_theErrorMatrixStore_Pull", "gaus", -theGaussianPullFitRange, theGaussianPullFitRange);
0700
0701     for (unsigned int i = 0; i != 5; ++i) {
0702       for (unsigned int j = i; j < 5; ++j) {
0703         unsigned int iH = theErrorMatrixStore_Pull->index(i, j);
0704         TProfile3D* ij = theErrorMatrixStore_Pull->get(i, j);
0705         TProfile3D* ii = theErrorMatrixStore_Pull->get(i, i);
0706         TProfile3D* jj = theErrorMatrixStore_Pull->get(j, j);
0707         if (!ij) {
0708           edm::LogError(theCategory) << i << " " << j << " not valid indexes. TProfile3D not found.";
0709           continue;
0710         }
0711         if (!ii) {
0712           edm::LogError(theCategory) << i << " " << i << " not valid indexes. TProfile3D not found.";
0713           continue;
0714         }
0715         if (!jj) {
0716           edm::LogError(theCategory) << j << " " << j << " not valid indexes. TProfile3D not found.";
0717           continue;
0718         }
0719
0720         unsigned int pTBin = (ij->GetNbinsX());
0721         unsigned int etaBin = (ij->GetNbinsY());
0722         unsigned int phiBin = (ij->GetNbinsZ());
0723
0724         //analyze the histograms now
0725         for (unsigned int iPt = 0; iPt < pTBin; ++iPt) {
0726           for (unsigned int iEta = 0; iEta < etaBin; iEta++) {
0727             for (unsigned int iPhi = 0; iPhi < phiBin; iPhi++) {
0728               double pt = ij->GetXaxis()->GetBinCenter(iPt + 1);
0729               double eta = ij->GetYaxis()->GetBinCenter(iEta + 1);
0730               double phi = ij->GetZaxis()->GetBinCenter(iPhi + 1);
0731
0732               TH1ptr& theH = (theHist_array_pull[iH])[index(ij, iPt, iEta, iPhi)];
0733               LogDebug(theCategory) << "extracting for: " << pt << " " << eta << " " << phi
0734                                     << "at index: " << index(ij, iPt, iEta, iPhi)
0735                                     << "while maxIndex is: " << maxIndex(ij) << "\n named: " << theH->GetName();
0736               double value = 0;
0737               if (i != j) {
0738                 //off diag term. not implemented yet. fill with ones
0739                 ij->Fill(pt, eta, phi, 1.);
0740               } else {
0741                 //diag term. perform a gaussian core fit (-2,2) to determine value
0742                 ((TH1*)theH)->Fit(f->GetName(), "qnL", "");
0743                 value = f->GetParameter(2);
0744                 LogDebug(theCategory) << "scale factor is: " << value;
0745                 ii->Fill(pt, eta, phi, value);
0746               }
0747
0748               //please free the memory !!!
0749               LogDebug(theCategory) << "freing memory of: " << theH->GetName();
0750               theH->Delete();
0751               theH = 0;
0752             }
0753           }
0754         }  //end of loop over the pt,eta,phi
0755       }
0756     }
0757
0758     theErrorMatrixStore_Pull->close();
0759   }
0760
0761   //close the file with all the plots in it.
0762   if (thePlotFile) {
0763     thePlotFile->Close();
0764   }
0765 }