Reference Manual

Go to the documentation of this file.
00001 #include "RecoParticleFlow/PFClusterTools/interface/TreeUtility.h"
00002 #include "TBranch.h"
00003 #include "TTree.h"
00004 #include "DataFormats/ParticleFlowReco/interface/Calibratable.h"
00005 #include <cmath>
00006 #include <fstream>
00007 #include <iostream>
00008 #include <TF1.h>
00009 
00010 using namespace pftools;
00011 TreeUtility::TreeUtility() {
00012 }
00013 
00014 TreeUtility::~TreeUtility() {
00015 }
00016 
00017 double deltaR(double eta1, double eta2, double phi1, double phi2) {
00018         return sqrt(pow(eta1 - eta2, 2) + pow(phi1 - phi2, 2));
00019 }
00020 
00021 unsigned TreeUtility::getCalibratablesFromRootFile(TChain& tree,
00022                 std::vector<Calibratable>& toBeFilled) {
00023 
00024         //      f.cd("extraction");
00025         //      TTree* tree = (TTree*) f.Get("extraction/Extraction");
00026         //      if (tree == 0) {
00027         //              PFToolsException me("Couldn't open tree!");
00028         //              throw me;
00029         //      }
00030         //      std::cout << "Successfully opened file. Getting branches..."<< std::endl;
00031         CalibratablePtr calib_ptr(new Calibratable());
00032         //TBranch* calibBr = tree.GetBranch("Calibratable");
00033         //spwBr->SetAddress(&spw);
00034         tree.SetBranchAddress("Calibratable", &calib_ptr);
00035 
00036         std::cout << "Looping over tree's "<< tree.GetEntries() << " entries...\n";
00037         for (unsigned entries(0); entries < tree.GetEntries(); entries++) {
00038                 tree.GetEntry(entries);
00039                 Calibratable c(*calib_ptr);
00040                 if (c.cands_num_ == 1 && (c.cluster_ecal_.size() +   c.cluster_hcal_.size()) > 0)
00041                         toBeFilled.push_back(c);
00042         }
00043         std::cout << "Done." << std::endl;
00044         return tree.GetEntries();
00045 
00046 }
00047 
00048 void TreeUtility::dumpCaloDataToCSV(TChain& tree, std::string csvFilename, double range, bool gaus) {
00049 
00050         CalibratablePtr calib_ptr(new Calibratable());
00051 
00052         tree.SetBranchAddress("Calibratable", &calib_ptr);
00053         ofstream csvFile;
00054         csvFile.open(csvFilename.c_str());
00055 
00056         std::cout << "Looping over tree's "<< tree.GetEntries() << " entries...\n";
00057         unsigned writes(0);
00058         TFile freq("freq.root", "recreate");
00059         TH1F frequencies("f", "f", 50, 0, range);
00060         TF1 g("g", "gaus(0)");
00061         g.FixParameter(1, range/2.0);
00062         g.FixParameter(0, 1),
00063         g.FixParameter(2, range/4.0);
00064 
00065         for (unsigned entries(0); entries < tree.GetEntries(); entries++) {
00066                 tree.GetEntry(entries);
00067                 Calibratable c(*calib_ptr);
00068                 bool veto(false);
00069 
00070                 //Check vetos as usual
00071                 if (c.cands_num_ > 1)
00072                         veto = true;
00073                  if(c.cluster_ecal_.size() == 0  && c.cluster_hcal_.size() == 0)
00074                         veto = true;
00075                  if(!veto) {
00076                         if(frequencies.GetBinContent(static_cast<int>(floor(c.sim_energyEvent_)) + 1) < (3000) * g.Eval(c.sim_energyEvent_) ) {
00077                                 frequencies.Fill(static_cast<int>(floor(c.sim_energyEvent_)));
00078                                 c.recompute();
00079                                 csvFile << c.sim_energyEvent_ << "\t";
00080                                 /*
00081                                 csvFile << c.sim_energyEcal_ << "\t";
00082                                 csvFile << c.sim_energyHcal_ << "\t";
00083 
00084 
00085                                 csvFile << c.cluster_energyEcal_/range << "\t";
00086                                 csvFile << c.cluster_energyHcal_/range << "\t";
00087 
00088                                 CaloWindow newEcalWindow(c.cluster_meanEcal_.eta_, c.cluster_meanEcal_.phi_, 5, 0.01, 3);
00089                                 const std::vector<CalibratableElement>& ecal = c.cluster_ecal_;
00090                                 std::vector<CalibratableElement>::const_iterator cit = ecal.begin();
00091                                 for(; cit != ecal.end(); ++cit) {
00092                                         const CalibratableElement& hit = *cit;
00093                                         bool added = newEcalWindow.addHit(hit.eta_, hit.phi_, hit.energy_);
00094                                         if(!added)
00095                                                 veto = true;
00096                                 }
00097                                 */
00098 
00099                                 csvFile << fabs(c.sim_eta_/2) << "\n";
00100 
00101                                 ++writes;
00102                         }
00103                 }
00104 
00105 
00106         }
00107         frequencies.Print("frequencies.eps");
00108         frequencies.Write();
00109         freq.Close();
00110         std::cout << "Closing file " << csvFilename << " with " << writes << " entries.\n" << std::endl;
00111 
00112         csvFile.close();
00113 
00114 }
00115 
00116 unsigned TreeUtility::getParticleDepositsDirectly(TChain& sourceChain,
00117                 std::vector<ParticleDepositPtr>& toBeFilled, CalibrationTarget target,
00118                 DetectorElementPtr offset, DetectorElementPtr ecal,
00119                 DetectorElementPtr hcal, bool includeOffset) {
00120 
00121         CalibratablePtr calib_ptr(new Calibratable());
00122         sourceChain.SetBranchAddress("Calibratable", &calib_ptr);
00123         std::cout << __PRETTY_FUNCTION__ << std::endl;
00124         std::cout << "WARNING: Using fabs() for eta value assignments!\n";
00125         std::cout << "Cutting on > 1 PFCandidate.\n";
00126         std::cout << "Looping over tree's "<< sourceChain.GetEntries() << " entries...\n";
00127         //neither of these two are supported yet
00128         if (target == UNDEFINED || target == PFELEMENT)
00129                 return 0;
00130         unsigned count(0);
00131 
00132         for (unsigned entries(0); entries < sourceChain.GetEntries(); entries++) {
00133                 sourceChain.GetEntry(entries);
00134                 Calibratable c(*calib_ptr);
00135 
00136                 ParticleDepositPtr pd(new ParticleDeposit());
00137                 bool veto(false);
00138                 if (c.sim_isMC_) {
00139                         pd->setTruthEnergy(c.sim_energyEvent_);
00140                         pd->setEta(fabs(c.sim_eta_));
00141                         pd->setPhi(c.sim_phi_);
00142                         //TODO:: sort this out
00143                         if (c.sim_energyEvent_== 0)
00144                                 veto = true;
00145                 }
00146                 if (c.tb_isTB_) {
00147                         pd->setTruthEnergy(c.sim_energyEvent_);
00148                         pd->setEta(c.tb_eta_);
00149                         pd->setPhi(c.tb_phi_);
00150                         veto = false;
00151                 }
00152 
00153                 if (c.cands_num_ > 1)
00154                         veto = true;
00155 
00156                 std::cout << "WARNING: HARD CUT ON 100 GeV SIM PARTICLES!\n";
00157                 if(c.sim_energyEvent_ > 100)
00158                         veto = true;
00159 
00160                 if (target == CLUSTER) {
00161                         if (c.cluster_ecal_.size() == 0  && c.cluster_hcal_.size() ==0)
00162                                 veto = true;
00163                         //                      if (c.cluster_numEcal_ > 1|| c.cluster_numHcal_ > 1)
00164                         //                              veto = true;
00165                         //TODO: using fabs for eta! WARNING!!!
00166                         Deposition decal(ecal, fabs(c.cluster_meanEcal_.eta_),
00167                                         c.cluster_meanEcal_.phi_, c.cluster_energyEcal_, 0);
00168                         Deposition dhcal(hcal, fabs(c.cluster_meanHcal_.eta_),
00169                                         c.cluster_meanHcal_.phi_, c.cluster_energyHcal_, 0);
00170                         Deposition doffset(offset, fabs(c.cluster_meanEcal_.eta_),
00171                                         c.cluster_meanEcal_.phi_, 0.001, 0);
00172 
00173                         pd->addTruthDeposition(decal);
00174                         pd->addRecDeposition(decal);
00175 
00176                         pd->addTruthDeposition(dhcal);
00177                         pd->addRecDeposition(dhcal);
00178 
00179                         if (includeOffset) {
00180                                 pd->addTruthDeposition(doffset);
00181                                 pd->addRecDeposition(doffset);
00182                         }
00183 
00184                 }
00185 
00186                 else if (target == PFCANDIDATE) {
00187                         //                      if(c.cands_num_ != 1)
00188                         //                              veto = true;
00189                         Deposition decal(ecal, c.cand_eta_, c.cand_phi_,
00190                                         c.cand_energyEcal_, 0);
00191                         Deposition dhcal(hcal, c.cand_eta_, c.cand_phi_,
00192                                         c.cand_energyHcal_, 0);
00193                         Deposition doffset(offset, c.cand_eta_, c.cand_phi_, 1.0, 0);
00194 
00195                         pd->addTruthDeposition(decal);
00196                         pd->addTruthDeposition(dhcal);
00197                         pd->addRecDeposition(decal);
00198                         pd->addRecDeposition(dhcal);
00199 
00200                         if (includeOffset) {
00201                                 pd->addTruthDeposition(doffset);
00202                                 pd->addRecDeposition(doffset);
00203                         }
00204                 }
00205 
00206                 else if (target == RECHIT) {
00207                         if (c.rechits_ecal_.size() == 0&& c.rechits_hcal_.size() == 0)
00208                                 veto = true;
00209                         Deposition decal(ecal, c.rechits_meanEcal_.eta_,
00210                                         c.rechits_meanEcal_.phi_, c.rechits_meanEcal_.energy_
00211                                                         * c.rechits_ecal_.size(), 0);
00212                         Deposition dhcal(hcal, c.rechits_meanHcal_.eta_,
00213                                         c.rechits_meanHcal_.phi_, c.rechits_meanHcal_.energy_
00214                                                         * c.rechits_hcal_.size(), 0);
00215                         Deposition doffset(offset, c.rechits_meanEcal_.eta_,
00216                                         c.rechits_meanEcal_.phi_, 1.0, 0);
00217 
00218                         pd->addTruthDeposition(decal);
00219                         pd->addTruthDeposition(dhcal);
00220                         pd->addRecDeposition(decal);
00221                         pd->addRecDeposition(dhcal);
00222 
00223                         if (includeOffset) {
00224                                 pd->addTruthDeposition(doffset);
00225                                 pd->addRecDeposition(doffset);
00226                         }
00227 
00228                 }
00229                 if (!veto)
00230                         toBeFilled.push_back(pd);
00231 
00232                 ++count;
00233         }
00234 
00235         return toBeFilled.size();
00236 }
00237 
00238 unsigned TreeUtility::convertCalibratablesToParticleDeposits(
00239                 const std::vector<Calibratable>& input,
00240                 std::vector<ParticleDepositPtr>& toBeFilled, CalibrationTarget target,
00241                 DetectorElementPtr offset, DetectorElementPtr ecal,
00242                 DetectorElementPtr hcal, bool includeOffset) {
00243 
00244         std::cout << __PRETTY_FUNCTION__ << std::endl;
00245         std::cout << "WARNING: Using fabs() for eta value assignments!\n";
00246         std::cout << "Input Calibratable has size "<< input.size() << "\n";
00247         std::cout << "Cutting on > 1 PFCandidate.\n";
00248 
00249         //neither of these two are supported yet
00250         if (target == UNDEFINED || target == PFELEMENT)
00251                 return 0;
00252         unsigned count(0);
00253         for (std::vector<Calibratable>::const_iterator cit = input.begin(); cit
00254                         != input.end(); ++cit) {
00255                 Calibratable c = *cit;
00256                 ParticleDepositPtr pd(new ParticleDeposit());
00257                 bool veto(false);
00258                 if (c.sim_isMC_) {
00259                         pd->setTruthEnergy(c.sim_energyEvent_);
00260                         pd->setEta(fabs(c.sim_eta_));
00261                         pd->setPhi(c.sim_phi_);
00262                         //TODO:: sort this out
00263                         if (c.sim_energyEvent_== 0)
00264                                 veto = true;
00265                 }
00266                 if (c.tb_isTB_) {
00267                         pd->setTruthEnergy(c.sim_energyEvent_);
00268                         pd->setEta(c.tb_eta_);
00269                         pd->setPhi(c.tb_phi_);
00270                         veto = false;
00271                 }
00272 
00273                 if (c.cands_num_ > 1)
00274                         veto = true;
00275 
00276                 if (target == CLUSTER) {
00277                         if (c.cluster_ecal_.size() == 0&& c.cluster_hcal_.size() ==0)
00278                                 veto = true;
00279                         //                      if (c.cluster_numEcal_ > 1|| c.cluster_numHcal_ > 1)
00280                         //                              veto = true;
00281                         //TODO: using fabs for eta! WARNING!!!
00282                         Deposition decal(ecal, fabs(c.cluster_meanEcal_.eta_),
00283                                         c.cluster_meanEcal_.phi_, c.cluster_energyEcal_, 0);
00284                         Deposition dhcal(hcal, fabs(c.cluster_meanHcal_.eta_),
00285                                         c.cluster_meanHcal_.phi_, c.cluster_energyHcal_, 0);
00286                         Deposition doffset(offset, fabs(c.cluster_meanEcal_.eta_),
00287                                         c.cluster_meanEcal_.phi_, 0.001, 0);
00288 
00289                         pd->addTruthDeposition(decal);
00290                         pd->addRecDeposition(decal);
00291 
00292                         pd->addTruthDeposition(dhcal);
00293                         pd->addRecDeposition(dhcal);
00294 
00295                         if (includeOffset) {
00296                                 pd->addTruthDeposition(doffset);
00297                                 pd->addRecDeposition(doffset);
00298                         }
00299 
00300                 }
00301 
00302                 else if (target == PFCANDIDATE) {
00303                         //                      if(c.cands_num_ != 1)
00304                         //                              veto = true;
00305                         Deposition decal(ecal, c.cand_eta_, c.cand_phi_,
00306                                         c.cand_energyEcal_, 0);
00307                         Deposition dhcal(hcal, c.cand_eta_, c.cand_phi_,
00308                                         c.cand_energyHcal_, 0);
00309                         Deposition doffset(offset, c.cand_eta_, c.cand_phi_, 1.0, 0);
00310 
00311                         pd->addTruthDeposition(decal);
00312                         pd->addTruthDeposition(dhcal);
00313                         pd->addRecDeposition(decal);
00314                         pd->addRecDeposition(dhcal);
00315 
00316                         if (includeOffset) {
00317                                 pd->addTruthDeposition(doffset);
00318                                 pd->addRecDeposition(doffset);
00319                         }
00320                 }
00321 
00322                 else if (target == RECHIT) {
00323                         if (c.rechits_ecal_.size() == 0&& c.rechits_hcal_.size() == 0)
00324                                 veto = true;
00325                         Deposition decal(ecal, c.rechits_meanEcal_.eta_,
00326                                         c.rechits_meanEcal_.phi_, c.rechits_meanEcal_.energy_
00327                                                         * c.rechits_ecal_.size(), 0);
00328                         Deposition dhcal(hcal, c.rechits_meanHcal_.eta_,
00329                                         c.rechits_meanHcal_.phi_, c.rechits_meanHcal_.energy_
00330                                                         * c.rechits_hcal_.size(), 0);
00331                         Deposition doffset(offset, c.rechits_meanEcal_.eta_,
00332                                         c.rechits_meanEcal_.phi_, 1.0, 0);
00333 
00334                         pd->addTruthDeposition(decal);
00335                         pd->addTruthDeposition(dhcal);
00336                         pd->addRecDeposition(decal);
00337                         pd->addRecDeposition(dhcal);
00338 
00339                         if (includeOffset) {
00340                                 pd->addTruthDeposition(doffset);
00341                                 pd->addRecDeposition(doffset);
00342                         }
00343 
00344                 }
00345                 if (!veto)
00346                         toBeFilled.push_back(pd);
00347 
00348                 ++count;
00349         }
00350 
00351         return toBeFilled.size();
00352 
00353 }
00354
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_5/src/RecoParticleFlow/PFClusterTools/src/TreeUtility.cc
