File indexing completed on 2024-04-06 12:21:37
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
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
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
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
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
0126 for(std::vector<event_data>::const_iterator i = data_.begin(); i != data_.end(); ++i)
0127 {
0128
0129
0130
0131
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
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
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
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
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 }