test/Calib_RootScripts/L1RCTCalibrator.C

0001 #define L1RCTCalibrator_cxx
0002 #include "L1RCTCalibrator.h"
0003 #include <TH2.h>
0004 #include <TStyle.h>
0005 #include <TCanvas.h>
0006 #include <iomanip>
0007
0008 #include <fstream>
0009
0010 void L1RCTCalibrator::Loop()
0011 {
0012 //   In a ROOT session, you can do:
0013 //      Root > .L L1RCTCalibrator.C
0014 //      Root > L1RCTCalibrator t
0015 //      Root > t.GetEntry(12); // Fill t data members with entry number 12
0016 //      Root > t.Show();       // Show values of entry 12
0017 //      Root > t.Show(16);     // Read and show values of entry 16
0018 //      Root > t.Loop();       // Loop on all entries
0019 //
0020
0021 //     This is the loop skeleton where:
0022 //    jentry is the global entry number in the chain
0023 //    ientry is the entry number in the current Tree
0024 //  Note that the argument to GetEntry must be:
0025 //    jentry for TChain::GetEntry
0026 //    ientry for TTree::GetEntry and TBranch::GetEntry
0027 //
0028 //       To read only selected branches, Insert statements like:
0029 // METHOD1:
0030 //    fChain->SetBranchStatus("*",0);  // disable all branches
0031 //    fChain->SetBranchStatus("branchname",1);  // activate branchname
0032 // METHOD2: replace line
0033 //    fChain->GetEntry(jentry);       //read all branches
0034 //by  b_branchname->GetEntry(ientry); //read only this branch
0035    if (fChain == 0) return;
0036
0037    Long64_t nentries = fChain->GetEntriesFast();
0038
0039    Long64_t nbytes = 0, nb = 0;
0040    for (Long64_t jentry=0; jentry<nentries;jentry++) {
0041       Long64_t ientry = LoadTree(jentry);
0042       if (ientry < 0) break;
0043       nb = fChain->GetEntry(jentry);   nbytes += nb;
0044    }
0045 }
0046
0047 void L1RCTCalibrator::makeCalibration()
0048 {
0049   if (fChain == 0) return;
0050
0051    Long64_t nentries = fChain->GetEntries();
0052    unsigned prev = 0;
0053
0054    std::cout << nentries << " to Process!" << std::endl;
0055
0056    Long64_t nbytes = 0, nb = 0;
0057    for (Long64_t jentry=0; jentry<nentries;jentry++) {
0058       Long64_t ientry = LoadTree(jentry);
0059       if (ientry < 0) break;
0060
0061       int check = (int)(100*(((double)jentry + 1)/ nentries));
0062
0063       if(debug_ < 1 && ( jentry == 0 || check - prev == 1 ))
0064     {
0065       prev = check;
0066       std::cout << std::fixed << '\r' <<  std::setprecision(1) << 100*(((double)jentry + 1)/ nentries) << "% data caching finished!" << std::flush;
0067     }
0068
0069
0070       nb = fChain->GetEntry(jentry);
0071       nbytes += nb;
0072
0073       event_data temp;
0074
0075       temp.event = Event_event;
0076       temp.run = Event_run;
0077
0078       //std::cout << temp.event << ' ' << temp.run << ' ' << Generator_nGen << ' ' <<  Region_nRegions << ' ' << CaloTPG_nTPG << std::endl;
0079
0080       for(unsigned i = 0; i < Generator_nGen; ++i)
0081     {
0082       generator tgen;
0083       tgen.particle_type = Generator_particle_type[i];
0084       tgen.et = Generator_et[i];
0085       tgen.phi = Generator_phi[i];
0086       tgen.eta = Generator_eta[i];
0087       temp.gen_particles.push_back(tgen);
0088     }
0089       for(unsigned i = 0; i < Region_nRegions; ++i)
0090     {
0091       region treg;
0092       treg.linear_et = Region_linear_et[i];
0093       treg.ieta = Region_ieta[i];
0094       treg.iphi = Region_iphi[i];
0095       treg.eta = Region_eta[i];
0096       treg.phi = Region_phi[i];
0097       temp.regions.push_back(treg);
0098     }
0099       for(unsigned i = 0; i < CaloTPG_nTPG; ++i)
0100     {
0101       tpg ttpg;
0102       ttpg.ieta = CaloTPG_ieta[i];
0103       ttpg.iphi = CaloTPG_iphi[i];
0104       ttpg.eta = CaloTPG_eta[i];
0105       ttpg.phi = CaloTPG_phi[i];
0106       ttpg.ecalEt = CaloTPG_ecalEt[i];
0107       ttpg.hcalEt = CaloTPG_hcalEt[i];
0108       ttpg.ecalE = CaloTPG_ecalE[i];
0109       ttpg.hcalE = CaloTPG_hcalE[i];
0110       temp.tpgs.push_back(ttpg);
0111     }
0112       data_.push_back(temp);
0113    }
0114
0115    std::cout << "\nDONE CACHING DATA TO MEMORY" << std::endl;
0116
0117    //now that we have all the data stored in memory, let's process it.. quickly...
0118
0119    if(debug_ > -1) std::cout << "------------------postProcessing()-------------------\n";
0120    int iph[28] = {0}, ipi[28] = {0}, ipi2[28] = {0},  nEvents = 0, nUseful = 0;
0121
0122    prev = 0;
0123
0124    std::cout << "Phase One: Collect Underpants" << std::endl;
0125    // first event data loop, calibrate ecal, hcal with NI pions
0126    for(std::vector<event_data>::const_iterator i = data_.begin(); i != data_.end(); ++i)
0127      {
0128        /*
0129      for(int n = 0; n < 28; ++n)
0130      {
0131      std::cout << "Trigger Tower " <<  n + 1 << ": " << iph[n] << ' ' << ipi[n] << ' ' << ipi2[n] << std::endl;
0132      }
0133        */
0134
0135        int check = (int)(100*(((double)nEvents + 1)/ nentries));
0136
0137        if(debug_ < 1 && ( nEvents == 0 || check - prev == 1 ))
0138      {
0139        prev = check;
0140        std::cout << std::fixed << '\r' << std::setprecision(1) << 100*(((double)nEvents + 1)/ nentries) << "% finished with polynomial fit!" << std::flush;
0141      }
0142
0143        nEvents++;
0144        hEvent->Fill(i->event);
0145        hRun->Fill(i->run);
0146
0147        bool used_flag = false;
0148
0149        std::vector<generator>::const_iterator gen = i->gen_particles.begin();
0150        std::vector<generator> ovls = overlaps(i->gen_particles);
0151
0152
0153       if(debug_ > 0) std::cout << "Finding overlapping selected particles in this event!" << std::endl;
0154       if(debug_ > 0)
0155     {
0156       std::cout << "=======Overlapping accepted gen particles in event " << nEvents << "" << std::endl;
0157       for(std::vector<generator>::const_iterator ov = ovls.begin(); ov != ovls.end(); ++ov)
0158         std::cout << '\t' <<  ov->particle_type << " " << ov->et << " " << ov->eta << " " << ov->phi << "" << std::endl;
0159       std::cout << "======================================================" << std::endl;
0160     }
0161       if(debug_ > 0) std::cout << "Done finding overlaps!" << std::endl;
0162
0163       for(; gen != i->gen_particles.end(); ++gen)
0164     {
0165       if(gen->particle_type != 22 && abs(gen->particle_type) != 211 || find<generator>(*gen, ovls) != -1) continue;
0166       double regionsum = sumEt(gen->eta,gen->phi,i->regions);
0167       if(regionsum > 0.0)
0168         {
0169           std::vector<tpg> matchedTpgs = tpgsNear(gen->eta,gen->phi,i->tpgs);
0170           std::pair<double,double> matchedCentroid = std::make_pair(avgEta(matchedTpgs),avgPhi(matchedTpgs));
0171           std::pair<double,double> etAndDeltaR95 = showerSize(matchedTpgs);
0172
0173           if(gen->particle_type == 22)
0174         {
0175           if(!used_flag) used_flag = true;
0176
0177           std::pair<double,double> ecalEtandDeltaR95 = showerSize(matchedTpgs, .95, .5, true, false);
0178
0179           int sumieta;
0180           if(ecalEtandDeltaR95.first > 0)
0181             {
0182               etaBin(fabs(matchedCentroid.first), sumieta);
0183
0184               hPhotonDeltaEOverE_uncor[sumieta - 1]->Fill((gen->et - ecalEtandDeltaR95.first)/gen->et);
0185
0186               if(debug_ > 0)
0187             {
0188               int n_towers = tpgsNear(matchedCentroid.first, matchedCentroid.second, matchedTpgs, ecalEtandDeltaR95.second).size();
0189               std::cout << "TPG sum near gen Photon with nearby non-zero RCT region found:" << std::endl
0190                     << "Number of Towers  : " << n_towers << " "
0191                     << "\tCentroid Eta      : " << matchedCentroid.first << " "
0192                     << "\tCentroid Phi      : " << matchedCentroid.second << " "
0193                     << "\tDelta R 95  (h+e) : " << etAndDeltaR95.second << " "
0194                     << "\nDelta R 95  (e)   : " << ecalEtandDeltaR95.second << " "
0195                     << "\tTotal Et in Cone  : " << etAndDeltaR95.first << " "
0196                     << "\tEcal Et in Cone   : " << ecalEtandDeltaR95.first << " ";
0197             }
0198
0199               hPhotonDeltaR95[sumieta - 1]->Fill(etAndDeltaR95.second);
0200               gPhotonEtvsGenEt[sumieta - 1]->SetPoint(iph[sumieta - 1]++, ecalEtandDeltaR95.first, gen->et);
0201             }
0202         }
0203
0204           if(abs(gen->particle_type) == 211)
0205         {
0206           if(!used_flag) used_flag = true;
0207
0208           double ecal = sumEt(matchedCentroid.first, matchedCentroid.second, matchedTpgs, etAndDeltaR95.second, true, false);
0209           double hcal = sumEt(matchedCentroid.first, matchedCentroid.second, matchedTpgs, etAndDeltaR95.second, false, true);
0210
0211           int sumieta;
0212           etaBin(fabs(matchedCentroid.first), sumieta);
0213
0214           hPionDeltaEOverE_uncor[sumieta - 1]->Fill((gen->et - (ecal + hcal))/gen->et);
0215
0216           if( ecal < 1.0  && etAndDeltaR95.first > 0)
0217             {
0218               if(debug_ > 0)
0219             {
0220               int n_towers = tpgsNear(matchedCentroid.first, matchedCentroid.second, matchedTpgs, etAndDeltaR95.second).size();
0221              std::cout << "TPG sum near gen Charged Pion with nearby non-zero RCT region and little ECAL energy found:" << std::endl
0222                    << "Number of Towers  : " << n_towers << " "
0223                    << "\tCentroid Eta      : " << matchedCentroid.first << " "
0224                    << "\tCentroid Phi      : " << matchedCentroid.second << " "
0225                    << "\tDelta R 95  (h+e) : " << etAndDeltaR95.second << " "
0226                    << "\nTotal Et in Cone  : " << etAndDeltaR95.first << " "
0227                    << "\tEcal Et in Cone   : " << ecal << " "
0228                    << "\tHcal Et in Cone   : " << hcal << " ";
0229             }
0230
0231               hNIPionDeltaR95[sumieta - 1]->Fill(etAndDeltaR95.second);
0232               gNIPionEtvsGenEt[sumieta - 1]->SetPoint(ipi[sumieta - 1]++, hcal, gen->et);
0233             }
0234           else if(etAndDeltaR95.first > 0)
0235             {
0236               if(debug_ > 0)
0237             {
0238               int n_towers = tpgsNear(matchedCentroid.first, matchedCentroid.second, matchedTpgs, etAndDeltaR95.second).size();
0239               std::cout << "TPG sum near gen Charged Pion with nearby non-zero RCT region found:" << std::endl
0240                     << "Number of Towers  : " << n_towers << " "
0241                     << "\tCentroid Eta      : " << matchedCentroid.first << " "
0242                     << "\tCentroid Phi      : " << matchedCentroid.second << " "
0243                     << "\tDelta R 95  (h+e) : " << etAndDeltaR95.second << " "
0244                     << "\nTotal Et in Cone  : " << etAndDeltaR95.first << " "
0245                     << "\tEcal Et in Cone   : " << ecal << " "
0246                     << "\tHcal Et in Cone   : " << hcal << " ";
0247             }
0248
0249               hPionDeltaR95[sumieta - 1]->Fill(etAndDeltaR95.second);
0250               gPionEcalEtvsHcalEtvsGenEt[sumieta - 1]->SetPoint(ipi2[sumieta - 1]++, ecal, hcal, gen->et);
0251             }
0252         }
0253
0254           hGenPhivsTpgSumPhi->Fill(gen->phi, matchedCentroid.second);
0255           hGenEtavsTpgSumEta->Fill(gen->eta, matchedCentroid.first);
0256           hTpgSumEt->Fill(etAndDeltaR95.first);
0257           hTpgSumEta->Fill(matchedCentroid.first);
0258           hTpgSumPhi->Fill(matchedCentroid.second);
0259         }
0260     }
0261
0262       if(used_flag) ++nUseful;
0263     }
0264
0265   int pitot = 0, phtot = 0, pitot2 = 0;
0266   TF1 *photon[28] = {NULL};
0267   TF1 *NIpion_low[28] = {NULL};
0268   TF1 *NIpion_high[28] = {NULL};
0269
0270   TF2* eh_surf[28] = {NULL};
0271
0272   std::cout << "\nBegin Fit" << std::endl;
0273   for(int i = 0; i < 28; ++i)
0274     {
0275       std::cout << "Photon and NI/I Pion Counts: " << iph[i] << ' ' << ipi[i] << ' ' << ipi2[i] << std::endl;
0276       pitot += ipi[i];
0277       phtot += iph[i];
0278       pitot2 += ipi2[i];
0279
0280       if(iph[i] > 100)
0281     {
0282       std::cout << "Fitting Photons!" << std::endl;
0283       photon[i] = new TF1((TString("ecal_fit")+=i).Data(),"x**3++x**2++x",0,100);
0284       gPhotonEtvsGenEt[i]->Fit(photon[i],fitOpts().c_str());
0285
0286       ecal_[i][0] = photon[i]->GetParameter(0);
0287       ecal_[i][1] = photon[i]->GetParameter(1);
0288       ecal_[i][2] = photon[i]->GetParameter(2);
0289     }
0290       else
0291     {
0292       ecal_[i][0] = 0;
0293       ecal_[i][1] = 0;
0294       ecal_[i][2] = 1;
0295     }
0296       if(ipi[i] > 100)
0297     {
0298       std::cout << "Fitting NI Pions!" << std::endl;
0299       NIpion_low[i] = new TF1((TString("hcal_fit_low")+=i).Data(),"x**3++x**2++x",0,100);
0300       NIpion_high[i] = new TF1((TString("hcal_fit_high")+=i).Data(),"x**3++x**2++x",0,100);
0301       gNIPionEtvsGenEt[i]->Fit(NIpion_low[i],fitOpts().c_str());
0302       gNIPionEtvsGenEt[i]->Fit(NIpion_high[i],fitOpts().c_str());
0303
0304       hcal_[i][0] = NIpion_low[i]->GetParameter(0);
0305       hcal_[i][1] = NIpion_low[i]->GetParameter(1);
0306       hcal_[i][2] = NIpion_low[i]->GetParameter(2);
0307
0308       hcal_high_[i][0] = NIpion_high[i]->GetParameter(0);
0309       hcal_high_[i][1] = NIpion_high[i]->GetParameter(1);
0310       hcal_high_[i][2] = NIpion_high[i]->GetParameter(2);
0311     }
0312       else
0313     {
0314       hcal_[i][0] = hcal_[i][1] = hcal_high_[i][0] = hcal_high_[i][1] = 0;
0315       hcal_[i][2] = hcal_high_[i][2] = 1;
0316     }
0317       if(false && ipi2[i] >100 && NIpion_low[i] && NIpion_high[i] && photon[i])
0318     {
0319       std::cout << "Fitting Pions" << std::endl;
0320       eh_surf[i] = new TF2((TString("eh_surf")+=i).Data(),"x**3++x**2++x++y**3++y**2++y++x**2*y++y**2*x++x*y++x**3*y++y**3*x++x**2*y**2",0,23,0,23);
0321       eh_surf[i]->FixParameter(0,photon[i]->GetParameter(0));
0322       eh_surf[i]->FixParameter(1,photon[i]->GetParameter(1));
0323       eh_surf[i]->FixParameter(2,photon[i]->GetParameter(2));
0324       eh_surf[i]->FixParameter(3,NIpion_low[i]->GetParameter(0));
0325       eh_surf[i]->FixParameter(4,NIpion_low[i]->GetParameter(1));
0326       eh_surf[i]->FixParameter(5,NIpion_low[i]->GetParameter(2));
0327       gPionEcalEtvsHcalEtvsGenEt[i]->Fit(eh_surf[i],fitOpts().c_str());
0328
0329       cross_[i][0] = eh_surf[i]->GetParameter(6);
0330       cross_[i][1] = eh_surf[i]->GetParameter(7);
0331       cross_[i][2] = eh_surf[i]->GetParameter(8);
0332       cross_[i][3] = eh_surf[i]->GetParameter(9);
0333       cross_[i][4] = eh_surf[i]->GetParameter(10);
0334       cross_[i][5] = eh_surf[i]->GetParameter(11);
0335     }
0336       else
0337     {
0338       cross_[i][0] = cross_[i][1] = cross_[i][2] = cross_[i][3] = cross_[i][4] = cross_[i][5] = 0;
0339     }
0340     }
0341   std::cout << "End Fit" << std::endl;
0342
0343   unsigned e = 0;
0344
0345   std::cout << "Phase Two: ????" << std::endl;
0346   for(std::vector<event_data>::const_iterator i = data_.begin(); i != data_.end(); ++i)
0347     {
0348
0349       int check = (int)(100*(((double)e + 1)/ nentries));
0350
0351        if(debug_ < 1 && ( e == 0 || check - prev == 1 ))
0352      {
0353        prev = check;
0354        std::cout << std::fixed << '\r' << std::setprecision(1) << 100*(((double)e + 1)/ nentries) << "% smearing factors!" << std::flush;
0355      }
0356        ++e;
0357
0358       std::vector<generator>::const_iterator gen = i->gen_particles.begin();
0359
0360       std::vector<generator> ovls = overlaps(i->gen_particles);
0361
0362       if(debug_ > 0) std::cout << "Finding overlapping selected particles in this event!" << std::endl;
0363       if(debug_ > 0)
0364     {
0365       std::cout << "=======Overlapping accepted gen particles in event " << nEvents << "" << std::endl;
0366       for(std::vector<generator>::const_iterator ov = ovls.begin(); ov != ovls.end(); ++ov)
0367         std::cout << '\t' <<  ov->particle_type << " " << ov->et << " " << ov->eta << " " << ov->phi << "" << std::endl;
0368       std::cout << "======================================================" << std::endl;
0369     }
0370       if(debug_ > 0) std::cout << "Done finding overlaps!" << std::endl;
0371
0372       for(; gen != i->gen_particles.end(); ++gen)
0373     {
0374       if(gen->particle_type != 22 && abs(gen->particle_type) != 211 || find<generator>(*gen, ovls) != -1) continue;
0375
0376       double regionsum = sumEt(gen->eta,gen->phi,i->regions);
0377
0378       if(regionsum > 0.0)
0379         {
0380           std::vector<tpg> matchedTpgs = tpgsNear(gen->eta,gen->phi,i->tpgs);
0381           std::pair<double,double> matchedCentroid = std::make_pair(avgEta(matchedTpgs),avgPhi(matchedTpgs));
0382           std::pair<double,double> etAndDeltaR95 = showerSize(matchedTpgs);
0383
0384           if(gen->particle_type == 22)
0385         {
0386           std::pair<double,double> ecalEtandDeltaR95 = showerSize(matchedTpgs, .95, .5, true, false);
0387
0388           int sumieta;
0389
0390           etaBin(fabs(matchedCentroid.first), sumieta);
0391
0392           double ecal_c = sumEt(matchedCentroid.first, matchedCentroid.second, matchedTpgs, etAndDeltaR95.second, true, false, true);
0393
0394           double deltaeovere = (gen->et - ecal_c)/gen->et;
0395           hPhotonDeltaEOverE[sumieta - 1]->Fill(deltaeovere);
0396
0397         }
0398
0399           if(abs(gen->particle_type) == 211)
0400         {
0401           double et_c = sumEt(matchedCentroid.first, matchedCentroid.second, matchedTpgs, etAndDeltaR95.second, true, true, true);
0402
0403           int sumieta;
0404           etaBin(fabs(matchedCentroid.first), sumieta);
0405
0406           double deltaeovere = (gen->et - et_c)/gen->et;
0407           hPionDeltaEOverE[sumieta - 1]->Fill(deltaeovere);
0408         }
0409         }
0410     }
0411     }
0412
0413   std::cout << std::endl;
0414
0415   TF1* phGaus[28];
0416   TF1* piGaus[28];
0417
0418   for(int i = 0; i < 28; ++i)
0419     {
0420       double ph_peak  = hPhotonDeltaEOverE[i]->GetBinCenter(hPhotonDeltaEOverE[i]->GetMaximumBin());
0421       double ph_upper = ph_peak + hPhotonDeltaEOverE[i]->GetRMS();
0422       double ph_lower = ph_peak - hPhotonDeltaEOverE[i]->GetRMS();
0423
0424       double pi_peak  = hPionDeltaEOverE[i]->GetBinCenter(hPionDeltaEOverE[i]->GetMaximumBin());
0425       double pi_upper = pi_peak + hPionDeltaEOverE[i]->GetRMS();
0426       double pi_lower = pi_peak - hPionDeltaEOverE[i]->GetRMS();
0427
0428       phGaus[i] = new TF1(TString("phGaus")+=i,"gaus",ph_lower,ph_upper);
0429       piGaus[i] = new TF1(TString("piGaus")+=i,"gaus",pi_lower,pi_upper);
0430
0431       if(hPhotonDeltaEOverE[i]->GetEntries() > 100)
0432     {
0433       hPhotonDeltaEOverE[i]->Fit(phGaus[i],fitOpts().c_str());
0434       he_low_smear_[i] = 1/(1 - phGaus[i]->GetParameter(1));
0435     }
0436       else
0437     he_low_smear_[i] = 1;
0438
0439       if(hPionDeltaEOverE[i]->GetEntries() > 100)
0440     {
0441       hPionDeltaEOverE[i]->Fit(piGaus[i],fitOpts().c_str());
0442       he_high_smear_[i] = 1/(1 - piGaus[i]->GetParameter(1));
0443     }
0444       else
0445     he_high_smear_[i] = 1;
0446     }
0447
0448   for(int i = 0; i < 28; ++i)
0449     {
0450       if(phGaus[i]) delete phGaus[i];
0451       if(phGaus[i]) delete piGaus[i];
0452       if(photon[i]) delete photon[i];
0453       if(NIpion_low[i]) delete NIpion_low[i];
0454       if(NIpion_high[i]) delete NIpion_high[i];
0455       if(eh_surf[i]) delete eh_surf[i];
0456     }
0457
0458   e = 0;
0459
0460   std::cout << "Phase 3: Profit!" << std::endl;
0461   for(std::vector<event_data>::const_iterator i = data_.begin(); i != data_.end(); ++i)
0462     {
0463       int check = (int)(100*(((double)e + 1)/ nentries));
0464
0465        if(debug_ < 1 && ( e == 0 || check - prev == 1 ))
0466      {
0467        prev = check;
0468        std::cout << std::fixed << '\r' << std::setprecision(1) << 100*(((double)nEvents + 1)/ nentries) << "% finished with diagnostics!" << std::flush;
0469      }
0470        ++e;
0471
0472       std::vector<generator>::const_iterator gen = i->gen_particles.begin();
0473
0474       std::vector<generator> ovls = overlaps(i->gen_particles);
0475
0476       if(debug_ > 0) std::cout << "Finding overlapping selected particles in this event!" << std::endl;
0477       if(debug_ > 0)
0478     {
0479       std::cout << "=======Overlapping accepted gen particles in event " << nEvents << "" << std::endl;
0480       for(std::vector<generator>::const_iterator ov = ovls.begin(); ov != ovls.end(); ++ov)
0481         std::cout << '\t' <<  ov->particle_type << " " << ov->et << " " << ov->eta << " " << ov->phi << "" << std::endl;
0482       std::cout << "======================================================" << std::endl;
0483     }
0484       if(debug_ > 0) std::cout << "Done finding overlaps!" << std::endl;
0485
0486       for(; gen != i->gen_particles.end(); ++gen)
0487     {
0488       if(gen->particle_type != 22 && abs(gen->particle_type) != 211 || find<generator>(*gen, ovls) != -1) continue;
0489
0490       double regionsum = sumEt(gen->eta,gen->phi,i->regions);
0491
0492       if(regionsum > 0.0)
0493         {
0494           std::vector<tpg> matchedTpgs = tpgsNear(gen->eta,gen->phi,i->tpgs);
0495           std::pair<double,double> matchedCentroid = std::make_pair(avgEta(matchedTpgs),avgPhi(matchedTpgs));
0496           std::pair<double,double> etAndDeltaR95 = showerSize(matchedTpgs);
0497
0498           double ecal = sumEt(matchedCentroid.first, matchedCentroid.second, matchedTpgs, etAndDeltaR95.second, true, false);
0499           double hcal = sumEt(matchedCentroid.first, matchedCentroid.second, matchedTpgs, etAndDeltaR95.second, false, true);
0500           double et_c = sumEt(matchedCentroid.first, matchedCentroid.second, matchedTpgs, etAndDeltaR95.second, true, true, true);
0501
0502           int sumieta;
0503           etaBin(fabs(matchedCentroid.first), sumieta);
0504
0505           if(hcal/(ecal + hcal) > 0.05 && etAndDeltaR95.first > 0)
0506         et_c *= ( (he_high_smear_[sumieta - 1] != -999) ? he_high_smear_[sumieta - 1] : 1 );
0507           else
0508         et_c *= ( (he_low_smear_[sumieta - 1] != -999) ? he_low_smear_[sumieta - 1] : 1 );
0509
0510           double deltaeovere = (gen->et - et_c)/gen->et;
0511
0512           if(gen->particle_type == 22)
0513         hPhotonDeltaEOverE_cor[sumieta - 1]->Fill(deltaeovere);
0514
0515           if(abs(gen->particle_type) == 211)
0516         {
0517           ecal = sumEt(matchedCentroid.first, matchedCentroid.second, matchedTpgs, etAndDeltaR95.second, true, false);
0518           hcal = sumEt(matchedCentroid.first, matchedCentroid.second, matchedTpgs, etAndDeltaR95.second, false, true, true);
0519           hPionDeltaEOverE_cor[sumieta - 1]->Fill(deltaeovere);
0520         }
0521         }
0522     }
0523     }
0524
0525   std::cout << std::endl;
0526
0527   // Now fill quick diagnostic plots.
0528   for(int i = 0; i < 25; ++i)
0529     {
0530       Double_t upper = hPhotonDeltaEOverE_cor[i]->GetBinCenter(hPhotonDeltaEOverE_cor[i]->GetMaximumBin()) + hPhotonDeltaEOverE_cor[i]->GetRMS();
0531       Double_t lower = hPhotonDeltaEOverE_cor[i]->GetBinCenter(hPhotonDeltaEOverE_cor[i]->GetMaximumBin()) - hPhotonDeltaEOverE_cor[i]->GetRMS();
0532
0533       hPhotonDeltaEOverE_cor[i]->Fit("gaus",
0534                      "QELM",
0535                      "",
0536                      lower,
0537                      upper);
0538
0539       hPhotonDeltaEtPeakvsEtaBinAllEt_c->SetBinContent(i+1,hPhotonDeltaEOverE_cor[i]->GetFunction("gaus")->GetParameter(1));
0540       hPhotonDeltaEtPeakvsEtaBinAllEt_c->SetBinError(i+1,hPhotonDeltaEOverE_cor[i]->GetFunction("gaus")->GetParError(1));
0541
0542       upper = hPhotonDeltaEOverE_uncor[i]->GetBinCenter(hPhotonDeltaEOverE_uncor[i]->GetMaximumBin()) + hPhotonDeltaEOverE_uncor[i]->GetRMS();
0543       lower = hPhotonDeltaEOverE_uncor[i]->GetBinCenter(hPhotonDeltaEOverE_uncor[i]->GetMaximumBin()) - hPhotonDeltaEOverE_uncor[i]->GetRMS();
0544
0545       hPhotonDeltaEOverE_uncor[i]->Fit("gaus",
0546                  "QELM",
0547                  "",
0548                  lower,
0549                  upper);
0550
0551       hPhotonDeltaEtPeakvsEtaBinAllEt_uc->SetBinContent(i+1,hPhotonDeltaEOverE_uncor[i]->GetFunction("gaus")->GetParameter(1));
0552       hPhotonDeltaEtPeakvsEtaBinAllEt_uc->SetBinError(i+1,hPhotonDeltaEOverE_uncor[i]->GetFunction("gaus")->GetParError(1));
0553
0554       upper = hPionDeltaEOverE_cor[i]->GetBinCenter(hPionDeltaEOverE_cor[i]->GetMaximumBin()) + hPionDeltaEOverE_cor[i]->GetRMS();
0555       lower = hPionDeltaEOverE_cor[i]->GetBinCenter(hPionDeltaEOverE_cor[i]->GetMaximumBin()) - hPionDeltaEOverE_cor[i]->GetRMS();
0556
0557       hPionDeltaEOverE_cor[i]->Fit("gaus",
0558                    "QELM",
0559                    "",
0560                    lower,
0561                    upper);
0562
0563       hPionDeltaEtPeakvsEtaBinAllEt_c->SetBinContent(i+1,hPionDeltaEOverE_cor[i]->GetFunction("gaus")->GetParameter(1));
0564       hPionDeltaEtPeakvsEtaBinAllEt_c->SetBinError(i+1,hPionDeltaEOverE_cor[i]->GetFunction("gaus")->GetParError(1));
0565
0566       upper = hPionDeltaEOverE_uncor[i]->GetBinCenter(hPionDeltaEOverE_uncor[i]->GetMaximumBin()) + hPionDeltaEOverE_uncor[i]->GetRMS();
0567       lower = hPionDeltaEOverE_uncor[i]->GetBinCenter(hPionDeltaEOverE_uncor[i]->GetMaximumBin()) - hPionDeltaEOverE_uncor[i]->GetRMS();
0568
0569       hPionDeltaEOverE_uncor[i]->Fit("gaus",
0570                      "QELM",
0571                      "",
0572                      lower,
0573                      upper);
0574
0575       hPionDeltaEtPeakvsEtaBinAllEt_uc->SetBinContent(i+1,hPionDeltaEOverE_uncor[i]->GetFunction("gaus")->GetParameter(1));
0576       hPionDeltaEtPeakvsEtaBinAllEt_uc->SetBinError(i+1,hPionDeltaEOverE_uncor[i]->GetFunction("gaus")->GetParError(1));
0577     }
0578
0579   if(debug_ > -1)
0580     std::cout << "Completed L1RCTGenCalibrator::postProcessing() with the following results:" << std::endl
0581           << "Total Events       : " << nEvents << "" << std::endl
0582           << "Useful Events      : " << nUseful << "" << std::endl
0583           << "Number of NI Pions : " << pitot << "" << std::endl
0584           << "Number of I Pions  : " << pitot2 << "" << std::endl
0585           << "Number of Photons  : " << phtot << "" << std::endl;
0586
0587   ofstream out;
0588   out.open("RCTCalibration_cff.py");
0589   WriteHistos();
0590   printCfFragment(out);
0591   out.close();
0592 }
0593
0594 void L1RCTCalibrator::BookHistos()
0595 {
0596   double deltaRbins[29];
0597
0598   for(int i = 0; i < 28; ++i)
0599     {
0600       double delta_r, eta, phi;
0601       phiValue(0,phi);
0602       etaValue(i+1,eta);
0603       deltaR(0,0,eta,phi,delta_r);
0604
0605       deltaRbins[i+1] = delta_r;
0606     }
0607
0608   putHist(hEvent = new TH1F("hEvent","Event Number",10000,0,10000));
0609   putHist(hRun = new TH1F("hRun","Run Number",10000,0,10000));
0610
0611   putHist(hGenPhi = new TH1F("hGenPhi","Generator #phi",72, 0, 2*M_PI));
0612   putHist(hGenEta = new TH1F("hGenEta","Generator #eta",100,-6,6));
0613   putHist(hGenEt = new TH1F("hGenEt","Generator E_{T}", 1000,0,500));
0614   putHist(hGenEtSel = new TH1F("hGenEtSel","Generator E_{T} Selected for Calibration",1000,0,500));
0615
0616   putHist(hRCTRegionEt = new TH1F("hRCTRegionEt","RCT Region E_{T} used in Calibration", 1000, 0, 500));
0617   putHist(hRCTRegionPhi = new TH1F("hRCTRegionPhi","RCT Region #phi", 72, 0, 2*M_PI));
0618   putHist(hRCTRegionEta = new TH1F("hRCTRegionEta","RCT Region #eta", 23, -5, 5));
0619
0620   putHist(hTpgSumEt = new TH1F("hTpgSumEt","TPG Sum E_{T}",1000,0,500));
0621   putHist(hTpgSumEta = new TH1F("hTpgSumEta","TPG Sum #eta Centroid", 57, -5,5));
0622   putHist(hTpgSumPhi = new TH1F("hTpgSumPhi","TPG Sum #phi Centroid", 72, 0, 2*M_PI));
0623
0624   putHist(hGenPhivsTpgSumPhi = new TH2F("hGenPhivsTpgSumPhi","Generator Particle #phi vs. TPG Sum #phi Centroid",72,0,2*M_PI,72,0,2*M_PI));
0625   putHist(hGenEtavsTpgSumEta = new TH2F("hGenEtavsTpgSumEta","Generator Particle #eta vs. TPG Sum #eta Centroid",57,-5,5,57,-5,5));
0626   putHist(hGenPhivsRegionPhi = new TH2F("hGenPhivsRegionPhi","Generator Particle #phi vs. RCT Region #phi",72,0,2*M_PI,72,0,2*M_PI));
0627   putHist(hGenEtavsRegionEta = new TH2F("hGenEtavsRegionEta","Generator Particle #eta vs. RCT Region #eta",23,-5,5,23,-5,5));
0628
0629   for(int i = 0; i < 28; ++i)
0630     {
0631       putHist(hPhotonDeltaR95[i] = new TH1F(TString("hPhotonDeltaR95") += i, TString("#gamma #DeltaR Containing 95% of E_{T} in #eta Bin: ") +=i, 28, deltaRbins));
0632       putHist(hNIPionDeltaR95[i] = new TH1F(TString("hNIPionDeltaR95") += i, TString("NI #pi^{#pm} #DeltaR Containing 95% of E_{T} in #eta Bin: ") +=i, 28, deltaRbins));
0633       putHist(hPionDeltaR95[i] = new TH1F(TString("hPionDeltaR95") += i, TString("#pi^{#pm} #DeltaR Containing 95% of E_{T} in #eta Bin: ") +=i, 28, deltaRbins));
0634
0635
0636       putHist(gPhotonEtvsGenEt[i] = new TGraphAsymmErrors());
0637       gPhotonEtvsGenEt[i]->SetName(TString("gPhotonEtvsGenEt") += i);
0638       gPhotonEtvsGenEt[i]->SetTitle(TString("#gamma TPG Sum E_{T} vs. Generator E_{T}, #eta Bin: ") += i);
0639
0640       putHist(gNIPionEtvsGenEt[i] = new TGraphAsymmErrors());
0641       gNIPionEtvsGenEt[i]->SetName(TString("gNIPionEtvsGenEt") += i);
0642       gNIPionEtvsGenEt[i]->SetTitle(TString("#pi^{#pm} with no/small ECAL deposit TPG Sum E_{T} vs. Generator E_{T}, #eta Bin: ") += i);
0643
0644       putHist(gPionEcalEtvsHcalEtvsGenEt[i] = new TGraph2DErrors());
0645       gPionEcalEtvsHcalEtvsGenEt[i]->SetName(TString("gPionEcalEtvsHcalEtvsGenEt") += i);
0646       gPionEcalEtvsHcalEtvsGenEt[i]->SetTitle(TString("#pi^{#pm} Ecal E_{T} vs Hcal E_{T} vs Generator E_{T}, #eta Bin: ") += i);
0647
0648       putHist(hPhotonDeltaEOverE[i] = new TH1F(TString("gPhotonDeltaEOverE")+=i,TString("Polynomial-fit-corrected #gamma #frac{#DeltaE_{T}}{E_{T}} in #eta Bin: ")+=i,
0649                            199,-2,2));
0650
0651       putHist(hPionDeltaEOverE[i] = new TH1F(TString("gPionDeltaEOverE")+=i,TString("Polynomial-fit-corrected #pi^{#pm} #frac{#DeltaE_{T}}{E_{T}} in #eta Bin: ")+=i,
0652                          199,-2,2));
0653
0654       putHist(hPhotonDeltaEOverE_cor[i] = new TH1F(TString("gPhotonDeltaEOverE_cor")+=i,TString("Corrected #gamma #frac{#DeltaE_{T}}{E_{T}} in #eta Bin: ")+=i,
0655                            199,-2,2));
0656       putHist(hPionDeltaEOverE_cor[i] = new TH1F(TString("gPionDeltaEOverE_cor")+=i,TString("Corrected #pi^{#pm} #frac{#DeltaE_{T}}{E_{T}} in #eta Bin: ")+=i,
0657                          199,-2,2));
0658
0659       putHist(hPhotonDeltaEOverE_uncor[i] = new TH1F(TString("gPhotonDeltaEOverE_uncor")+=i,TString("Uncorrected #gamma #frac{#DeltaE_{T}}{E_{T}} in #eta Bin: ")+=i,
0660                              199,-2,2));
0661       putHist(hPionDeltaEOverE_uncor[i] = new TH1F(TString("gPionDeltaEOverE_uncor")+=i,TString("Uncorrected #pi^{#pm} #frac{#DeltaE_{T}}{E_{T}} in #eta Bin: ")+=i,
0662                            199,-2,2));
0663
0664     }
0665
0666   for(int i = 0; i < 12; ++i)
0667     {
0668       putHist(hDeltaEtPeakvsEtaBin_uc[i] = new TH1F(TString("hDeltaEtPeakvsEtaBin_uc") += i,
0669                             TString("Uncorrected #DeltaE_{T} Peak vs. #eta Bin in E_{T} Bin ") += i,
0670                             25,0.5,25.5));
0671       putHist(hDeltaEtPeakvsEtaBin_c[i] = new TH1F(TString("hDeltaEtPeakvsEtaBin_c") += i,
0672                            TString("Corrected #DeltaE_{T} Peak vs. #eta Bin in E_{T} Bin ") += i,
0673                            25,0.5,25.5));
0674       putHist(hDeltaEtPeakRatiovsEtaBin[i] = new TH1F(TString("hEtPeakRatiovsEtaBin") += i,
0675                               TString("#frac{Corrected #E_{T}}{Uncorrected E_{T}} vs. #eta Bin in E_{T} Bin ") += i,
0676                               25,0.5,25.5));
0677     }
0678
0679   putHist(hPhotonDeltaEtPeakvsEtaBinAllEt_uc = new TH1F("hPhotonDeltaEtPeakvsEtaBinAllEt_uc",
0680                             "Photon Uncorrected #DeltaE_{T} vs. #eta Bin, All E_{T}",
0681                             25,0.5,25.5));
0682   putHist(hPhotonDeltaEtPeakvsEtaBinAllEt_c = new TH1F("hPhotonDeltaEtPeakvsEtaBinAllEt_c",
0683                                "Photon Corrected #DeltaE_{T} vs. #eta Bin, All E_{T}",
0684                                25,0.5,25.5));
0685   putHist(hPhotonDeltaEtPeakRatiovsEtaBinAllEt = new TH1F("hPhotonDeltaEtPeakRatiovsEtaBinAllEt",
0686                               "Photon #frac{Corrected E_{T}}{Uncorrected E_{T} vs. #eta Bin, All E_{T}",
0687                               25,0.5,25.5));
0688
0689   putHist(hPionDeltaEtPeakvsEtaBinAllEt_uc = new TH1F("hPionDeltaEtPeakvsEtaBinAllEt_uc",
0690                             "Pion Uncorrected #DeltaE_{T} vs. #eta Bin, All E_{T}",
0691                             25,0.5,25.5));
0692   putHist(hPionDeltaEtPeakvsEtaBinAllEt_c = new TH1F("hPionDeltaEtPeakvsEtaBinAllEt_c",
0693                                "Pion Corrected #DeltaE_{T} vs. #eta Bin, All E_{T}",
0694                                25,0.5,25.5));
0695   putHist(hPionDeltaEtPeakRatiovsEtaBinAllEt = new TH1F("hPionDeltaEtPeakRatiovsEtaBinAllEt",
0696                               "Pion #frac{Corrected E_{T}}{Uncorrected E_{T} vs. #eta Bin, All E_{T}",
0697                               25,0.5,25.5));
0698 }
0699
0700 void L1RCTCalibrator::WriteHistos()
0701 {
0702   output_->cd();
0703
0704   for(std::vector<TObject*>::const_iterator i = hists_.begin(); i != hists_.end(); ++i)
0705     (*i)->Write();
0706
0707   output_->Write();
0708   output_->Close();
0709 }
0710
0711
0712 double L1RCTCalibrator::uniPhi(const double& phi) const
0713 {
0714   double result = ((phi < 0) ? phi + 2*M_PI : phi);
0715   while(result > 2*M_PI) result -= 2*M_PI;
0716   return result;
0717 }
0718
0719 //calculate Delta R between two (eta,phi) coordinates
0720 void L1RCTCalibrator::deltaR(const double& eta1, double phi1,
0721                  const double& eta2, double phi2,double& dr) const
0722 {
0723   double deta2 = std::pow(eta1-eta2,2.);
0724   double dphi2 = std::pow(uniPhi(phi1)-uniPhi(phi2),2.);
0725
0726   dr = std::sqrt(deta2 + dphi2);
0727 }
0728
0729 void L1RCTCalibrator::etaBin(const double& veta, int& ieta) const
0730 {
0731   double absEta = fabs(veta);
0732
0733   if(absEta < maxEtaBarrel_)
0734     {
0735       ieta = static_cast<int>((absEta+0.000001)/deltaEtaBarrel_) + 1;
0736     }
0737   else
0738     {
0739       double temp = absEta - maxEtaBarrel_;
0740       int i = 0;
0741       while(temp > -0.0000001 && i < 8)
0742     {
0743       temp -= endcapEta_[i++];
0744     }
0745       ieta = 20 + i;
0746     }
0747   ieta = ((veta < 0) ? -ieta : ieta);
0748 }
0749
0750 void L1RCTCalibrator::etaValue(const int& ieta, double& veta) const
0751 {
0752   int absEta = abs(ieta);
0753
0754   if(absEta <= 20)
0755     {
0756       veta = (absEta-1)*0.0870 + 0.0870/2.;
0757     }
0758   else
0759     {
0760       int offset = abs(ieta) - 21;
0761       veta = maxEtaBarrel_;
0762       for(int i = 0; i < offset; ++i)
0763     veta += endcapEta_[i];
0764       veta += endcapEta_[offset]/2.;
0765     }
0766   veta = ((ieta < 0) ? -veta : veta);
0767 }
0768
0769 void L1RCTCalibrator::phiBin(double vphi, int& iphi) const
0770 {
0771   iphi = static_cast<int>(uniPhi(vphi)/deltaPhi_);
0772 }
0773
0774 void L1RCTCalibrator::phiValue(const int& iphi, double& vphi) const
0775 {
0776   vphi = iphi*deltaPhi_ + deltaPhi_/2.;
0777 }
0778
0779 bool L1RCTCalibrator::sanityCheck() const
0780 {
0781   for(int i = 1; i <= 28; ++i)
0782     {
0783       int j, l;
0784       double p, q;
0785       etaValue(i,p);
0786       etaBin(p,j);
0787       etaValue(-i,q);
0788       etaBin(q,l);
0789       if( i != j || -i != l)
0790     {
0791       if(debug_ > -1)
0792         std::cout << i <<  "\t" << p << "\t" << j << "\t"
0793               << -i << "\t" << q << "\t" << l << std::endl;
0794       return false;
0795     }
0796     }
0797   for(int i = 0; i < 72; ++i)
0798     {
0799       int j;
0800       double p;
0801       phiValue(i,p);
0802       phiBin(p,j);
0803       if(i != j)
0804     {
0805       if(debug_ > -1)
0806         std::cout << i << "\t" << p << "\t" << j << std::endl;
0807       return false;
0808     }
0809     }
0810   return true;
0811 }
0812
0813 // energy, deltaR
0814 std::pair<double,double> L1RCTCalibrator::showerSize(const std::vector<tpg>& tp, const double frac, const double& max_dr,
0815                              const bool& ecal, const bool& hcal) const
0816 {
0817   double c_eta = avgEta(tp), c_phi = avgPhi(tp);
0818   double result = 0.0,
0819     e_max = sumEt(c_eta, c_phi, tp, max_dr, ecal, hcal);
0820
0821   double dr_iter = 0.0;
0822
0823   do{
0824     result = sumEt(c_eta, c_phi, tp, dr_iter, ecal, hcal);
0825     dr_iter += 0.01;
0826   }while(result/e_max < frac);
0827
0828   return std::make_pair(result,dr_iter);
0829 }
0830
0831 double L1RCTCalibrator::sumEt(const double& eta, const double& phi, const std::vector<tpg>& tp, const double& dr,
0832                   const bool& ecal, const bool& hcal, const bool& c, const double& crossover) const
0833 {
0834   double delta_r, sum = 0.0;
0835
0836   for(std::vector<tpg>::const_iterator i = tp.begin(); i != tp.end(); ++i)
0837     {
0838       deltaR(eta,phi,i->eta,i->phi,delta_r);
0839
0840       if(delta_r < dr)
0841     {
0842       if(c)
0843         {
0844           int etabin = abs(i->ieta) - 1;
0845           if(i->ecalE > .5 && ecal)
0846         {
0847           if(ecal_[etabin][0] != -999 && ecal_[etabin][1] != -999 && ecal_[etabin][2] != -999)
0848             sum += (ecal_[etabin][0]*std::pow(i->ecalEt,3.) +
0849                 ecal_[etabin][1]*std::pow(i->ecalEt,2.) +
0850                 ecal_[etabin][2]*i->ecalEt);
0851           else
0852             sum += i->ecalEt;
0853         }
0854           if(i->hcalE > .5 && hcal)
0855         {
0856           double crossterm = 0.0, hcal_c = 0.0;
0857           if(i->ecalEt + i->hcalEt < crossover)
0858             {
0859               if(cross_[etabin][0] != -999 && cross_[etabin][1] != -999 && cross_[etabin][2] != -999 &&
0860              cross_[etabin][3] != -999 && cross_[etabin][4] != -999 && cross_[etabin][5] != -999 &&
0861              hcal_[etabin][0] != -999 && hcal_[etabin][1] != -999 && hcal_[etabin][2] != -999)
0862             {
0863               crossterm = (cross_[etabin][0]*std::pow(i->ecalEt,2)*i->hcalEt +
0864                        cross_[etabin][1]*std::pow(i->hcalEt,2)*i->ecalEt +
0865                        cross_[etabin][2]*i->ecalEt*i->hcalEt +
0866                        cross_[etabin][3]*std::pow(i->ecalEt,3)*i->hcalEt +
0867                        cross_[etabin][4]*std::pow(i->hcalEt,3)*i->ecalEt +
0868                        cross_[etabin][5]*std::pow(i->ecalEt,2)*std::pow(i->hcalEt,2));
0869               hcal_c = (hcal_[etabin][0]*std::pow(i->hcalEt,3.) +
0870                     hcal_[etabin][1]*std::pow(i->hcalEt,2.) +
0871                     hcal_[etabin][2]*i->hcalEt);
0872             }
0873               else
0874             hcal_c = i->hcalEt;
0875             }
0876           else
0877             {
0878               if(hcal_high_[etabin][0] != -999 && hcal_high_[etabin][1] != -999 && hcal_high_[etabin][2] != -999)
0879             {
0880               hcal_c = (hcal_high_[etabin][0]*std::pow(i->hcalEt,3.) +
0881                     hcal_high_[etabin][1]*std::pow(i->hcalEt,2.) +
0882                     hcal_high_[etabin][2]*i->hcalEt);
0883             }
0884               else
0885             hcal_c = i->hcalEt;
0886             }
0887           sum += hcal_c + crossterm;
0888         }
0889         }
0890       else
0891         {
0892           if(i->ecalE > .5 && ecal) sum += i->ecalEt;
0893           if(i->hcalE > .5 && hcal) sum += i->hcalEt;
0894         }
0895     }
0896     }
0897   return sum;
0898 }
0899
0900 double L1RCTCalibrator::sumEt(const double& eta, const double& phi, const std::vector<region>& regs, const double& dr) const
0901 {
0902   double sum = 0.0, delta_r;
0903
0904   for(std::vector<region>::const_iterator i = regs.begin(); i != regs.end(); ++i)
0905     {
0906       deltaR(eta,phi,i->eta,i->phi,delta_r);
0907
0908       if(delta_r < dr) sum += i->linear_et*.5;
0909     }
0910   return sum;
0911 }
0912
0913 double L1RCTCalibrator::avgPhi(const std::vector<tpg>& t) const
0914 {
0915   double n = 0.0, d = 0.0;
0916
0917   for(std::vector<tpg>::const_iterator i = t.begin(); i != t.end(); ++i)
0918     {
0919       n = (i->ecalEt + i->hcalEt)*i->phi;
0920       d = i->ecalEt + i->hcalEt;
0921     }
0922
0923   return n/d;
0924 }
0925
0926 double L1RCTCalibrator::avgEta(const std::vector<tpg>& t) const
0927 {
0928   double n = 0.0, d = 0.0;
0929
0930   for(std::vector<tpg>::const_iterator i = t.begin(); i != t.end(); ++i)
0931     {
0932       double temp_eta;
0933       etaValue(i->ieta, temp_eta);
0934       n = (i->ecalEt + i->hcalEt)*temp_eta;
0935       d = i->ecalEt + i->hcalEt;
0936     }
0937   return n/d;
0938 }
0939
0940 std::vector<tpg> L1RCTCalibrator::tpgsNear(const double& eta, const double& phi, const std::vector<tpg>& tpgs,
0941                        const double& dr) const
0942 {
0943   std::vector<tpg> result;
0944
0945   for(std::vector<tpg>::const_iterator i = tpgs.begin(); i != tpgs.end(); ++i)
0946     {
0947       double delta_r;
0948       deltaR(eta,phi,i->eta,i->phi,delta_r);
0949
0950       if(delta_r < dr) result.push_back(*i);
0951     }
0952   return result;
0953 }
0954
0955 //This prints out a nicely formatted .cfi file to be included in RCTConfigProducer.cfi
0956 void L1RCTCalibrator::printCfFragment(std::ostream& out) const
0957 {
0958   double* p = NULL;
0959
0960   out.flush();
0961
0962   out << ((python_) ? "import FWCore.ParameterSet.Config as cms\n\nrct_calibration = cms.PSet(" : "block rct_calibration = {") << std::endl;
0963   for(int i = 0; i < 6; ++i)
0964     {
0965       switch(i)
0966     {
0967     case 0:
0968       p = const_cast<double*>(reinterpret_cast<const double*>(ecal_));
0969       out << ((python_) ? "\tecal_calib_Lindsey = cms.vdouble(" : "\tvdouble ecal_calib_Lindsey = {") << std::endl;
0970       break;
0971     case 1:
0972       p = const_cast<double*>(reinterpret_cast<const double*>(hcal_));
0973       out << ((python_) ? "\thcal_calib_Lindsey = cms.vdouble(" : "\tvdouble hcal_calib_Lindsey = {") << std::endl;
0974       break;
0975     case 2:
0976       p = const_cast<double*>(reinterpret_cast<const double*>(hcal_high_));
0977       out << ((python_) ? "\thcal_high_calib_Lindsey = cms.vdouble(" : "\tvdouble hcal_high_calib_Lindsey = {") << std::endl;
0978       break;
0979     case 3:
0980       p = const_cast<double*>(reinterpret_cast<const double*>(cross_));
0981       out << ((python_) ? "\tcross_terms_Lindsey = cms.vdouble(" : "\tvdouble cross_terms_Lindsey = {") << std::endl;
0982       break;
0983     case 4:
0984       p = const_cast<double*>(reinterpret_cast<const double*>(he_low_smear_));
0985       out << ((python_) ? "\tHoverE_low_Lindsey = cms.vdouble(" : "\tvdouble HoverE_low_Lindsey = {") << std::endl;
0986       break;
0987     case 5:
0988       p = const_cast<double*>(reinterpret_cast<const double*>(he_high_smear_));
0989       out << ((python_) ? "\tHoverE_high_Lindsey = cms.vdouble(" : "\tvdouble HoverE_high_Lindsey = {") << std::endl;
0990     };
0991
0992       for(int j = 0; j < 28; ++j)
0993     {
0994       if( p == reinterpret_cast<const double*>(ecal_) || p == reinterpret_cast<const double*>(hcal_) ||
0995           p == reinterpret_cast<const double*>(hcal_high_) )
0996         {
0997           double *q = p + 3*j;
0998           if(q[0] != -999 && q[1] != -999 && q[2] != -999)
0999         {
1000           out << "\t\t" << q[0] << ", " << q[1] << ", " << q[2];
1001           out << ((j==27) ? "" : ",") << std::endl;
1002         }
1003         }
1004       else if( p == reinterpret_cast<const double*>(cross_) )
1005         {
1006           double *q = p + 6*j;
1007           if(q[0] != -999 && q[1] != -999 && q[2] != -999 &&
1008          q[3] != -999 && q[4] != -999 && q[5] != -999)
1009         {
1010           out << "\t\t" << q[0] << ", " << q[1] << ", " << q[2] << ", "
1011               << q[3] << ", " << q[4] << ", " << q[5];
1012           out << ((j==27) ? "" : ",") << std::endl;
1013         }
1014         }
1015       else
1016         {
1017           double *q = p;
1018           if(q[j] != -999)
1019         out << "\t\t" << q[j] << ((j==27) ? "" : ",") << std::endl;
1020         }
1021     }
1022       if(python_)
1023     {
1024       out << ((i != 5) ? "\t)," : "\t)") << std::endl;
1025     }
1026       else
1027     out << "\t}" << std::endl;
1028     }
1029   out << ((python_) ? ")" : "}") << std::endl;
1030 }
1031
1032 std::vector<generator>
1033 L1RCTCalibrator::overlaps(const std::vector<generator>& v) const
1034 {
1035   std::vector<generator> result;
1036
1037   if(v.begin() == v.end()) return result;
1038
1039   for(std::vector<generator>::const_iterator i = v.begin(); i != v.end() - 1; ++i)
1040     {
1041       for(std::vector<generator>::const_iterator j = i + 1; j != v.end(); ++j)
1042     {
1043       double delta_r;
1044
1045       deltaR(i->eta, i->phi, j->eta, j->phi, delta_r);
1046
1047       if(delta_r < .5)
1048         {
1049           if(find<generator>(*i, result) == -1)
1050         result.push_back(*i);
1051           if(find<generator>(*j, result) == -1)
1052         result.push_back(*j);
1053         }
1054     }
1055     }
1056
1057   return result;
1058 }