MomentumScaleCalibration/plugins/MuScleFitPlotter.cc

0001 //  \class MuScleFitPlotter
0002 //  Plotter for simulated,generated and reco info of muons
0003 //
0004 //  \author  C.Mariotti, S.Bolognesi - INFN Torino / T.Dorigo, M.De Mattia - INFN Padova
0005 // revised S. Casasso, E. Migliore - UniTo & INFN Torino
0006 //
0007 // ----------------------------------------------------------------------------------
0008
0009 #include "MuScleFitPlotter.h"
0010 #include "Histograms.h"
0011 #include "MuScleFitUtils.h"
0012 #include "FWCore/Framework/interface/MakerMacros.h"
0013 #include "FWCore/Framework/interface/Frameworkfwd.h"
0014 #include "FWCore/Framework/interface/Event.h"
0015 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0016 //#include "SimGeneral/HepPDTRecord/interface/ParticleDataTable.h"
0017 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0018 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
0019 #include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
0020 #include "DataFormats/TrackReco/interface/Track.h"
0021 #include "DataFormats/MuonReco/interface/Muon.h"
0022 #include "DataFormats/MuonReco/interface/MuonFwd.h"
0023 #include "DataFormats/Candidate/interface/LeafCandidate.h"
0024 #include "RecoMuon/TrackingTools/interface/MuonPatternRecoDumper.h"
0025 #include "RecoMuon/TrackingTools/interface/MuonServiceProxy.h"
0026 #include <CLHEP/Vector/LorentzVector.h>
0027
0028 #include "TFile.h"
0029 #include "TTree.h"
0030 #include "TMinuit.h"
0031 #include <vector>
0032
0033 // Constructor
0034 // ----------
0035 MuScleFitPlotter::MuScleFitPlotter(std::string theGenInfoRootFileName) {
0036   outputFile = new TFile(theGenInfoRootFileName.c_str(), "RECREATE");
0037   fillHistoMap();
0038 }
0039
0040 MuScleFitPlotter::~MuScleFitPlotter() {
0041   outputFile->cd();
0042   writeHistoMap();
0043   outputFile->Close();
0044 }
0045
0046 // Find and store in histograms the generated resonance and muons
0047 // --------------------------------------------------------------
0048 void MuScleFitPlotter::fillGen(const reco::GenParticleCollection& genParticles, bool PATmuons) {
0049   //  bool prova = false;
0050   //Loop on generated particles
0051   std::pair<reco::Particle::LorentzVector, reco::Particle::LorentzVector> muFromRes;
0052   reco::Particle::LorentzVector genRes;
0053
0054   int mothersFound[] = {0, 0, 0, 0, 0, 0};
0055
0056   for (reco::GenParticleCollection::const_iterator mcIter = genParticles.begin(); mcIter != genParticles.end();
0057        ++mcIter) {
0058     int status = mcIter->status();
0059     int pdgId = std::abs(mcIter->pdgId());
0060     //Check if it's a resonance
0061     if (status == 2 &&
0062         (pdgId == 23 || pdgId == 443 || pdgId == 100443 || pdgId == 553 || pdgId == 100553 || pdgId == 200553)) {
0063       genRes = mcIter->p4();
0064       // std::cout << "mother's mother = " << mcIter->mother()->pdgId() << std::endl;
0065       if (pdgId == 23)
0066         mapHisto["hGenResZ"]->Fill(genRes);
0067       else if (pdgId == 443 || pdgId == 100443)
0068         mapHisto["hGenResJPsi"]->Fill(genRes);
0069       else if (pdgId == 553 || pdgId == 100553 || pdgId == 200553)
0070         mapHisto["hGenResUpsilon1S"]->Fill(genRes);
0071     }
0072     //Check if it's a muon from a resonance
0073     if (status == 1 && pdgId == 13 && !PATmuons) {
0074       int momPdgId = std::abs(mcIter->mother()->pdgId());
0075       if (momPdgId == 23 || momPdgId == 443 || momPdgId == 100443 || momPdgId == 553 || momPdgId == 100553 ||
0076           momPdgId == 200553) {
0077         if (momPdgId == 23)
0078           mothersFound[0] = 1;
0079         if (momPdgId == 443 || momPdgId == 100443)
0080           mothersFound[5] = 1;
0081         if (momPdgId == 553 || momPdgId == 100553 || momPdgId == 200553)
0082           mothersFound[3] = 1;
0083         mapHisto["hGenMu"]->Fill(mcIter->p4());
0084         std::cout << "genmu " << mcIter->p4() << std::endl;
0085         if (mcIter->charge() > 0) {
0086           muFromRes.first = mcIter->p4();
0087           // prova = true;
0088         } else
0089           muFromRes.second = mcIter->p4();
0090       }
0091     }  //if PATmuons you don't have the info of the mother !!! Here I assume is a JPsi
0092     if (status == 1 && pdgId == 13 && PATmuons) {
0093       mothersFound[5] = 1;
0094       mapHisto["hGenMu"]->Fill(mcIter->p4());
0095       std::cout << "genmu " << mcIter->p4() << std::endl;
0096       if (mcIter->charge() > 0) {
0097         muFromRes.first = mcIter->p4();
0098         // prova = true;
0099       } else
0100         muFromRes.second = mcIter->p4();
0101     }
0102   }
0103   //   if(!prova)
0104   //     std::cout<<"hgenmumu not found"<<std::endl;
0105
0106   if (mothersFound[0] == 1) {
0107     mapHisto["hGenMuMuZ"]->Fill(muFromRes.first + muFromRes.second);
0108     mapHisto["hGenResVSMuZ"]->Fill(muFromRes.first, genRes, 1);
0109     mapHisto["hGenResVSMuZ"]->Fill(muFromRes.second, genRes, -1);
0110   }
0111   if (mothersFound[3] == 1) {
0112     mapHisto["hGenMuMuUpsilon1S"]->Fill(muFromRes.first + muFromRes.second);
0113     mapHisto["hGenResVSMuUpsilon1S"]->Fill(muFromRes.first, genRes, 1);
0114     mapHisto["hGenResVSMuUpsilon1S"]->Fill(muFromRes.second, genRes, -1);
0115   }
0116   if (mothersFound[5] == 1) {
0117     mapHisto["hGenMuMuJPsi"]->Fill(muFromRes.first + muFromRes.second);
0118     mapHisto["hGenResVSMuJPsi"]->Fill(muFromRes.first, genRes, 1);
0119     mapHisto["hGenResVSMuJPsi"]->Fill(muFromRes.second, genRes, -1);
0120   }
0121
0122   mapHisto["hGenResVsSelf"]->Fill(genRes, genRes, 1);
0123 }
0124
0125 // Find and store in histograms the generated resonance and muons
0126 // --------------------------------------------------------------
0127 void MuScleFitPlotter::fillGen(const edm::HepMCProduct& evtMC, bool sherpaFlag_) {
0128   //Loop on generated particles
0129   const HepMC::GenEvent* Evt = evtMC.GetEvent();
0130   std::pair<reco::Particle::LorentzVector, reco::Particle::LorentzVector> muFromRes;
0131   reco::Particle::LorentzVector genRes;
0132
0133   int mothersFound[] = {0, 0, 0, 0, 0, 0};
0134
0135   if (sherpaFlag_) {
0136     for (HepMC::GenEvent::particle_const_iterator part = Evt->particles_begin(); part != Evt->particles_end(); part++) {
0137       if (std::abs((*part)->pdg_id()) == 13 && (*part)->status() == 1) {  //looks for muon in the final state
0138         bool fromRes = false;
0139         for (HepMC::GenVertex::particle_iterator mother = (*part)->production_vertex()->particles_begin(
0140                  HepMC::ancestors);  //loops on the mother of the final state muons
0141              mother != (*part)->production_vertex()->particles_end(HepMC::ancestors);
0142              ++mother) {
0143           unsigned int motherPdgId = (*mother)->pdg_id();
0144           if (motherPdgId == 13 && (*mother)->status() == 3)
0145             fromRes = true;
0146         }
0147         if (fromRes) {
0148           if ((*part)->pdg_id() == 13) {
0149             muFromRes.first = (lorentzVector(
0150                 (*part)->momentum().px(), (*part)->momentum().py(), (*part)->momentum().pz(), (*part)->momentum().e()));
0151           } else if ((*part)->pdg_id() == -13) {
0152             muFromRes.second = (lorentzVector(
0153                 (*part)->momentum().px(), (*part)->momentum().py(), (*part)->momentum().pz(), (*part)->momentum().e()));
0154           }
0155         }
0156       }
0157     }
0158     mapHisto["hGenResZ"]->Fill(muFromRes.first + muFromRes.second);
0159   } else {
0160     for (HepMC::GenEvent::particle_const_iterator part = Evt->particles_begin(); part != Evt->particles_end(); part++) {
0161       int status = (*part)->status();
0162       int pdgId = std::abs((*part)->pdg_id());
0163       //std::cout<<"PDG ID "<< (*part)->pdg_id() <<"    status "<< (*part)->status()
0164       //<<"   pt "<<(*part)->momentum().perp()<< "     eta  "<<(*part)->momentum().eta()<<std::endl    ;
0165       //Check if it's a resonance
0166       // For sherpa the resonance is not saved. The muons from the resonance can be identified
0167       // by having as mother a muon of status 3.
0168
0169       if (pdgId == 13 && status == 1) {
0170         if (status == 2 &&
0171             (pdgId == 23 || pdgId == 443 || pdgId == 100443 || pdgId == 553 || pdgId == 100553 || pdgId == 200553)) {
0172           genRes = reco::Particle::LorentzVector(
0173               (*part)->momentum().px(), (*part)->momentum().py(), (*part)->momentum().pz(), (*part)->momentum().e());
0174
0175           if (pdgId == 23)
0176             mapHisto["hGenResZ"]->Fill(genRes);
0177           if (pdgId == 443)
0178             mapHisto["hGenResJPsi"]->Fill(genRes);
0179           if (pdgId == 553) {
0180             // std::cout << "genRes mass = " << CLHEP::HepLorentzVector(genRes.x(),genRes.y(),genRes.z(),genRes.t()).m() << std::endl;
0181             mapHisto["hGenResUpsilon1S"]->Fill(genRes);
0182           }
0183         }
0184
0185         //Check if it's a muon from a resonance
0186         if (pdgId == 13 && status == 1) {
0187           bool fromRes = false;
0188           for (HepMC::GenVertex::particle_iterator mother =
0189                    (*part)->production_vertex()->particles_begin(HepMC::ancestors);
0190                mother != (*part)->production_vertex()->particles_end(HepMC::ancestors);
0191                ++mother) {
0192             int motherPdgId = (*mother)->pdg_id();
0193             if (motherPdgId == 23 || motherPdgId == 443 || motherPdgId == 100443 || motherPdgId == 553 ||
0194                 motherPdgId == 100553 || motherPdgId == 200553) {
0195               fromRes = true;
0196               if (motherPdgId == 23)
0197                 mothersFound[0] = 1;
0198               if (motherPdgId == 443)
0199                 mothersFound[3] = 1;
0200               if (motherPdgId == 553)
0201                 mothersFound[5] = 1;
0202             }
0203           }
0204
0205           if (fromRes) {
0206             mapHisto["hGenMu"]->Fill(reco::Particle::LorentzVector(
0207                 (*part)->momentum().px(), (*part)->momentum().py(), (*part)->momentum().pz(), (*part)->momentum().e()));
0208             mapHisto["hGenMuVSEta"]->Fill(reco::Particle::LorentzVector(
0209                 (*part)->momentum().px(), (*part)->momentum().py(), (*part)->momentum().pz(), (*part)->momentum().e()));
0210             if ((*part)->pdg_id() == -13)
0211               muFromRes.first = (reco::Particle::LorentzVector((*part)->momentum().px(),
0212                                                                (*part)->momentum().py(),
0213                                                                (*part)->momentum().pz(),
0214                                                                (*part)->momentum().e()));
0215             else
0216               muFromRes.second = (reco::Particle::LorentzVector((*part)->momentum().px(),
0217                                                                 (*part)->momentum().py(),
0218                                                                 (*part)->momentum().pz(),
0219                                                                 (*part)->momentum().e()));
0220           }
0221         }
0222       }
0223     }
0224   }
0225   if (mothersFound[0] == 1) {
0226     mapHisto["hGenMuMuZ"]->Fill(muFromRes.first + muFromRes.second);
0227     mapHisto["hGenResVSMuZ"]->Fill(muFromRes.first, genRes, 1);
0228     mapHisto["hGenResVSMuZ"]->Fill(muFromRes.second, genRes, -1);
0229   }
0230   if (mothersFound[3] == 1) {
0231     mapHisto["hGenMuMuUpsilon1S"]->Fill(muFromRes.first + muFromRes.second);
0232     mapHisto["hGenResVSMuUpsilon1S"]->Fill(muFromRes.first, genRes, 1);
0233     mapHisto["hGenResVSMuUpsilon1S"]->Fill(muFromRes.second, genRes, -1);
0234   }
0235   if (mothersFound[5] == 1) {
0236     mapHisto["hGenMuMuJPsi"]->Fill(muFromRes.first + muFromRes.second);
0237     mapHisto["hGenResVSMuJPsi"]->Fill(muFromRes.first, genRes, 1);
0238     mapHisto["hGenResVSMuJPsi"]->Fill(muFromRes.second, genRes, -1);
0239   }
0240   mapHisto["hGenResVsSelf"]->Fill(genRes, genRes, 1);
0241 }
0242
0243 // Find and store in histograms the simulated resonance and muons
0244 // --------------------------------------------------------------
0245 void MuScleFitPlotter::fillSim(edm::Handle<edm::SimTrackContainer> simTracks) {
0246   std::vector<SimTrack> simMuons;
0247
0248   //Loop on simulated tracks
0249   for (edm::SimTrackContainer::const_iterator simTrack = simTracks->begin(); simTrack != simTracks->end(); ++simTrack) {
0250     // Select the muons from all the simulated tracks
0251     if (std::abs((*simTrack).type()) == 13) {
0252       simMuons.push_back(*simTrack);
0253       mapHisto["hSimMu"]->Fill((*simTrack).momentum());
0254     }
0255   }
0256   mapHisto["hSimMu"]->Fill(simMuons.size());
0257
0258   // Recombine all the possible Z from simulated muons
0259   if (simMuons.size() >= 2) {
0260     for (std::vector<SimTrack>::const_iterator imu = simMuons.begin(); imu != simMuons.end(); ++imu) {
0261       for (std::vector<SimTrack>::const_iterator imu2 = imu + 1; imu2 != simMuons.end(); ++imu2) {
0262         if (imu == imu2)
0263           continue;
0264
0265         // Try all the pairs with opposite charge
0266         if (((*imu).charge() * (*imu2).charge()) < 0) {
0267           reco::Particle::LorentzVector Z = (*imu).momentum() + (*imu2).momentum();
0268           mapHisto["hSimMuPMuM"]->Fill(Z);
0269         }
0270       }
0271     }
0272
0273     // Plots for the best possible simulated resonance
0274     std::pair<SimTrack, SimTrack> simMuFromBestRes = MuScleFitUtils::findBestSimuRes(simMuons);
0275     reco::Particle::LorentzVector bestSimZ = (simMuFromBestRes.first).momentum() + (simMuFromBestRes.second).momentum();
0276     mapHisto["hSimBestRes"]->Fill(bestSimZ);
0277     if (std::abs(simMuFromBestRes.first.momentum().eta()) < 2.5 &&
0278         std::abs(simMuFromBestRes.second.momentum().eta()) < 2.5 && simMuFromBestRes.first.momentum().pt() > 2.5 &&
0279         simMuFromBestRes.second.momentum().pt() > 2.5) {
0280       mapHisto["hSimBestResVSMu"]->Fill(
0281           simMuFromBestRes.first.momentum(), bestSimZ, int(simMuFromBestRes.first.charge()));
0282       mapHisto["hSimBestResVSMu"]->Fill(
0283           simMuFromBestRes.second.momentum(), bestSimZ, int(simMuFromBestRes.second.charge()));
0284     }
0285   }
0286 }
0287
0288 // Find and store in histograms the RIGHT simulated resonance and muons
0289 // --------------------------------------------------------------
0290 void MuScleFitPlotter::fillGenSim(edm::Handle<edm::HepMCProduct> evtMC, edm::Handle<edm::SimTrackContainer> simTracks) {
0291   std::pair<reco::Particle::LorentzVector, reco::Particle::LorentzVector> simMuFromRes =
0292       MuScleFitUtils::findSimMuFromRes(evtMC, simTracks);
0293   //Fill resonance info
0294   reco::Particle::LorentzVector rightSimRes = (simMuFromRes.first) + (simMuFromRes.second);
0295   mapHisto["hSimRightRes"]->Fill(rightSimRes);
0296   /*if ((std::abs(simMuFromRes.first.Eta())<2.5 && std::abs(simMuFromRes.second.Eta())<2.5)
0297     && simMuFromRes.first.Pt()>2.5 && simMuFromRes.second.Pt()>2.5) {
0298   }*/
0299 }
0300
0301 // // Find and store in histograms the reconstructed resonance and muons
0302 // // --------------------------------------------------------------
0303 // void MuScleFitPlotter::fillRec(std::vector<reco::LeafCandidate>& muons)
0304 // {
0305 //   for(std::vector<reco::LeafCandidate>::const_iterator mu1 = muons.begin(); mu1!=muons.end(); mu1++){
0306 //     mapHisto["hRecMu"]->Fill(mu1->p4());
0307 //     mapHisto["hRecMuVSEta"]->Fill(mu1->p4());
0308 //     for(std::vector<reco::LeafCandidate>::const_iterator mu2 = muons.begin(); mu2!=muons.end(); mu2++){
0309 //       if (mu1->charge()<0 || mu2->charge()>0)
0310 //  continue;
0311 //       reco::Particle::LorentzVector Res (mu1->p4()+mu2->p4());
0312 //       mapHisto["hRecMuPMuM"]->Fill(Res);
0313 //     }
0314 //   }
0315 //   mapHisto["hRecMu"]->Fill(muons.size());
0316 // }
0317
0318 /// Used when running on the root tree containing preselected muon pairs
0319 void MuScleFitPlotter::fillRec(std::vector<MuScleFitMuon>& muons) {
0320   for (std::vector<MuScleFitMuon>::const_iterator mu1 = muons.begin(); mu1 != muons.end(); mu1++) {
0321     mapHisto["hRecMu"]->Fill(mu1->p4());
0322     mapHisto["hRecMuVSEta"]->Fill(mu1->p4());
0323     for (std::vector<MuScleFitMuon>::const_iterator mu2 = muons.begin(); mu2 != muons.end(); mu2++) {
0324       if (mu1->charge() < 0 || mu2->charge() > 0)
0325         continue;
0326       reco::Particle::LorentzVector Res(mu1->p4() + mu2->p4());
0327       mapHisto["hRecMuPMuM"]->Fill(Res);
0328     }
0329   }
0330   mapHisto["hRecMu"]->Fill(muons.size());
0331 }
0332
0333 /// Used when running on the root tree containing preselected muon pairs
0334 void MuScleFitPlotter::fillTreeRec(
0335     const std::vector<std::pair<reco::Particle::LorentzVector, reco::Particle::LorentzVector> >& savedPairs) {
0336   std::vector<std::pair<reco::Particle::LorentzVector, reco::Particle::LorentzVector> >::const_iterator muonPair =
0337       savedPairs.begin();
0338   for (; muonPair != savedPairs.end(); ++muonPair) {
0339     mapHisto["hRecMu"]->Fill(muonPair->first);
0340     mapHisto["hRecMuVSEta"]->Fill(muonPair->first);
0341     mapHisto["hRecMu"]->Fill(muonPair->second);
0342     mapHisto["hRecMuVSEta"]->Fill(muonPair->second);
0343     reco::Particle::LorentzVector Res(muonPair->first + muonPair->second);
0344     mapHisto["hRecMuPMuM"]->Fill(Res);
0345     mapHisto["hRecMu"]->Fill(savedPairs.size());
0346   }
0347 }
0348
0349 /**
0350  * Used when running on the root tree and there is genInfo. <br>
0351  * ATTENTION: since we do not have any id information when reading from the root tree, we always
0352  * fill the Z histograms by default.
0353  */
0354 void MuScleFitPlotter::fillTreeGen(
0355     const std::vector<std::pair<reco::Particle::LorentzVector, reco::Particle::LorentzVector> >& genPairs) {
0356   std::vector<std::pair<reco::Particle::LorentzVector, reco::Particle::LorentzVector> >::const_iterator genPair =
0357       genPairs.begin();
0358   for (; genPair != genPairs.end(); ++genPair) {
0359     reco::Particle::LorentzVector genRes(genPair->first + genPair->second);
0360     mapHisto["hGenResZ"]->Fill(genRes);
0361     mapHisto["hGenMu"]->Fill(genPair->first);
0362     mapHisto["hGenMuVSEta"]->Fill(genPair->first);
0363     mapHisto["hGenMuMuZ"]->Fill(genRes);
0364     mapHisto["hGenResVSMuZ"]->Fill(genPair->first, genRes, 1);
0365     mapHisto["hGenResVSMuZ"]->Fill(genPair->second, genRes, -1);
0366     mapHisto["hGenMuMuUpsilon1S"]->Fill(genRes);
0367     mapHisto["hGenResVSMuUpsilon1S"]->Fill(genPair->first, genRes, 1);
0368     mapHisto["hGenResVSMuUpsilon1S"]->Fill(genPair->second, genRes, -1);
0369     mapHisto["hGenMuMuJPsi"]->Fill(genRes);
0370     mapHisto["hGenResVSMuJPsi"]->Fill(genPair->first, genRes, 1);
0371     mapHisto["hGenResVSMuJPsi"]->Fill(genPair->second, genRes, -1);
0372     mapHisto["hGenResVsSelf"]->Fill(genRes, genRes, 1);
0373   }
0374 }
0375
0376 // Histogram booking
0377 // -----------------
0378 void MuScleFitPlotter::fillHistoMap() {
0379   // Generated Z and muons
0380   // ---------------------
0381   mapHisto["hGenResJPsi"] = new HParticle("hGenResJPsi", 3., 3.1);
0382   mapHisto["hGenResUpsilon1S"] = new HParticle("hGenResUpsilon1S", 9., 11.);
0383   mapHisto["hGenResZ"] = new HParticle("hGenResZ", 60., 120.);
0384   mapHisto["hGenMu"] = new HParticle("hGenMu");
0385   mapHisto["hGenMuVSEta"] = new HPartVSEta("hGenMuVSEta");
0386
0387   mapHisto["hGenMuMuJPsi"] = new HParticle("hGenMuMuJPsi", 3., 3.1);
0388   mapHisto["hGenResVSMuJPsi"] = new HMassVSPart("hGenResVSMuJPsi", 3., 3.1);
0389   mapHisto["hGenMuMuUpsilon1S"] = new HParticle("hGenMuMuUpsilon1S", 9., 11.);
0390   mapHisto["hGenResVSMuUpsilon1S"] = new HMassVSPart("hGenResVSMuUpsilon1S", 9., 11.);
0391   mapHisto["hGenMuMuZ"] = new HParticle("hGenMuMuZ", 60., 120.);
0392   mapHisto["hGenResVSMuZ"] = new HMassVSPart("hGenResVSMuZ", 60., 120.);
0393
0394   mapHisto["hGenResVsSelf"] = new HMassVSPart("hGenResVsSelf");
0395
0396   // Simulated resonance and muons
0397   // -----------------------------
0398   mapHisto["hSimMu"] = new HParticle("hSimMu");
0399
0400   mapHisto["hSimMuPMuM"] = new HParticle("hSimMuPMuM");
0401
0402   mapHisto["hSimBestMu"] = new HParticle("hSimBestMu");
0403   mapHisto["hSimBestRes"] = new HParticle("hSimBestRes");
0404   mapHisto["hSimBestResVSMu"] = new HMassVSPart("hSimBestResVSMu");
0405
0406   mapHisto["hSimRightRes"] = new HParticle("hSimRightZ");
0407
0408   // Reconstructed resonance and muons
0409   // -----------------------------
0410   mapHisto["hRecMu"] = new HParticle("hRecMu");
0411   mapHisto["hRecMuVSEta"] = new HPartVSEta("hRecMuVSEta");
0412   mapHisto["hRecMuPMuM"] = new HParticle("hRecMuPMuM");
0413 }
0414
0415 // Histogram saving
0416 // -----------------
0417 void MuScleFitPlotter::writeHistoMap() {
0418   outputFile->cd();
0419   for (std::map<std::string, Histograms*>::const_iterator histo = mapHisto.begin(); histo != mapHisto.end(); histo++) {
0420     (*histo).second->Write();
0421   }
0422 }