73 EventsAfterCuts =
new TH1F(
"EventsAfterCuts",
"Events After Cuts", 30, 0, 30);
76 e9 =
new TH1F(
"e9",
"E9 energy", 300, 0., 150.);
77 e25 =
new TH1F(
"e25",
"E25 energy", 300, 0., 150.);
78 scE =
new TH1F(
"scE",
"SC energy", 300, 0., 150.);
79 trP =
new TH1F(
"trP",
"Trk momentum", 300, 0., 150.);
80 EoP =
new TH1F(
"EoP",
"EoP", 600, 0., 3.);
81 EoP_all =
new TH1F(
"EoP_all",
"EoP_all", 600, 0., 3.);
83 calibs =
new TH1F(
"calib",
"Calibration constants", 800, 0.5, 2.);
84 calibsEndCapMinus =
new TH1F(
"calibEndCapMinus",
"Calibration constants EE-", 800, 0.5, 2.);
85 calibsEndCapPlus =
new TH1F(
"calibEndCapPlus",
"Calibration constants EE+", 800, 0.5, 2.);
87 e25OverScE =
new TH1F(
"e25OverscE",
"E25 / SC energy", 400, 0., 2.);
88 E25oP =
new TH1F(
"E25oP",
"E25 / P", 750, 0., 1.5);
90 Map =
new TH2F(
"Map",
"Nb Events in Crystal", 173, -86, 86, 362, 0, 361);
91 e9Overe25 =
new TH1F(
"e9Overe25",
"E9 / E25", 400, 0., 2.);
92 Map3Dcalib =
new TH2F(
"3Dcalib",
"3Dcalib", 173, -86, 86, 362, 0, 361);
96 MapCor1 =
new TH2F(
"MapCor1",
"Correlation E25/Pcalo versus E25/Pin", 100, 0., 5., 100, 0., 5.);
97 MapCor2 =
new TH2F(
"MapCor2",
"Correlation E25/Pcalo versus E/P", 100, 0., 5., 100, 0., 5.);
98 MapCor3 =
new TH2F(
"MapCor3",
"Correlation E25/Pcalo versus Pout/Pin", 100, 0., 5., 100, 0., 5.);
99 MapCor4 =
new TH2F(
"MapCor4",
"Correlation E25/Pcalo versus E25/highestP", 100, 0., 5., 100, 0., 5.);
100 MapCor5 =
new TH2F(
"MapCor5",
"Correlation E25/Pcalo versus Pcalo/Pout", 100, 0., 5., 100, 0., 5.);
101 MapCor6 =
new TH2F(
"MapCor6",
"Correlation Pout/Pin versus E25/Pin", 100, 0., 5., 100, 0., 5.);
102 MapCor7 =
new TH2F(
"MapCor7",
"Correlation Pout/Pin versus Pcalo/Pout", 100, 0., 5., 100, 0., 5.);
103 MapCor8 =
new TH2F(
"MapCor8",
"Correlation E25/Pin versus Pcalo/Pout", 100, 0., 5., 100, 0., 5.);
104 MapCor9 =
new TH2F(
"MapCor9",
"Correlation E25/Pcalo versus Eseed/Pout", 100, 0., 5., 100, 0., 5.);
105 MapCor10 =
new TH2F(
"MapCor10",
"Correlation Eseed/Pout versus Pout/Pin", 100, 0., 5., 100, 0., 5.);
106 MapCor11 =
new TH2F(
"MapCor11",
"Correlation Eseed/Pout versus E25/Pin", 100, 0., 5., 100, 0., 5.);
109 E25oPvsEta =
new TH2F(
"E25oPvsEta",
"E/P vs Eta", 173, -86, 86, 600, 0.7, 1.3);
110 E25oPvsEtaEndCapMinus =
new TH2F(
"E25oPvsEtaEndCapMinus",
"E/P vs R EE-", 100, 0, 100, 600, 0.7, 1.3);
111 E25oPvsEtaEndCapPlus =
new TH2F(
"E25oPvsEtaEndCapPlus",
"E/P vs R EE+", 100, 0, 100, 600, 0.7, 1.3);
113 PinMinPout =
new TH1F(
"PinMinPout",
"(Pin - Pout)/Pin", 600, -2.0, 2.0);
115 calibinter =
new TH1F(
"calibinter",
"internal calibration constants", 800, 0.5, 2.);
116 PinOverPout =
new TH1F(
"PinOverPout",
"pinOverpout", 600, 0., 3.);
117 eSeedOverPout =
new TH1F(
"eSeedOverPout",
"eSeedOverpout ", 600, 0., 3.);
121 calibinterEndCapMinus =
new TH1F(
"calibinterEndCapMinus",
"internal calibration constants", 800, 0.5, 2.);
122 calibinterEndCapPlus =
new TH1F(
"calibinterEndCapPlus",
"internal calibration constants", 800, 0.5, 2.);
129 Error1 =
new TH1F(
"Error1",
"DeltaP/Pin", 800, -1.0, 1.0);
130 Error2 =
new TH1F(
"Error2",
"DeltaP/Pout", 800, -1.0, 1.0);
131 Error3 =
new TH1F(
"Error3",
"DeltaP/Pcalo", 800, -1.0, 1.0);
132 eSeedOverPout2 =
new TH1F(
"eSeedOverPout2",
"eSeedOverpout (No Supercluster)", 600, 0., 4.);
133 hadOverEm =
new TH1F(
"hadOverEm",
"Had/EM distribution", 600, -2., 2.);
136 Map3DcalibNoCuts =
new TH2F(
"3DcalibNoCuts",
"3Dcalib (Before Cuts)", 173, -86, 86, 362, 0, 361);
137 e9NoCuts =
new TH1F(
"e9NoCuts",
"E9 energy (Before Cuts)", 300, 0., 150.);
138 e25NoCuts =
new TH1F(
"e25NoCuts",
"E25 energy (Before Cuts)", 300, 0., 150.);
139 scENoCuts =
new TH1F(
"scENoCuts",
"SC energy (Before Cuts)", 300, 0., 150.);
140 trPNoCuts =
new TH1F(
"trPNoCuts",
"Trk momentum (Before Cuts)", 300, 0., 150.);
141 EoPNoCuts =
new TH1F(
"EoPNoCuts",
"EoP (Before Cuts)", 600, 0., 3.);
142 calibsNoCuts =
new TH1F(
"calibNoCuts",
"Calibration constants (Before Cuts)", 800, 0., 2.);
143 e25OverScENoCuts =
new TH1F(
"e25OverscENoCuts",
"E25 / SC energy (Before Cuts)", 400, 0., 2.);
144 E25oPNoCuts =
new TH1F(
"E25oPNoCuts",
"E25 / P (Before Cuts)", 750, 0., 1.5);
145 MapEndCapMinus =
new TH2F(
"MapEndCapMinus",
"Nb Events in Crystal (EndCap)", 100, 0, 100, 100, 0, 100);
146 MapEndCapPlus =
new TH2F(
"MapEndCapPlus",
"Nb Events in Crystal (EndCap)", 100, 0, 100, 100, 0, 100);
147 e9Overe25NoCuts =
new TH1F(
"e9Overe25NoCuts",
"E9 / E25 (Before Cuts)", 400, 0., 2.);
148 PinOverPoutNoCuts =
new TH1F(
"PinOverPoutNoCuts",
"pinOverpout (Before Cuts)", 600, 0., 3.);
149 eSeedOverPoutNoCuts =
new TH1F(
" eSeedOverPoutNoCuts",
"eSeedOverpout (Before Cuts) ", 600, 0., 4.);
150 PinMinPoutNoCuts =
new TH1F(
"PinMinPoutNoCuts",
"(Pin - Pout)/Pin (Before Cuts)", 600, -2.0, 2.0);
154 calibinterNoCuts =
new TH1F(
"calibinterNoCuts",
"internal calibration constants", 2000, 0.5, 2.);
157 new TH2F(
"MapCor1NoCuts",
"Correlation E25/PatCalo versus E25/Pin (Before Cuts)", 100, 0., 5., 100, 0., 5.);
159 new TH2F(
"MapCor2NoCuts",
"Correlation E25/PatCalo versus E/P (Before Cuts)", 100, 0., 5., 100, 0., 5.);
161 new TH2F(
"MapCor3NoCuts",
"Correlation E25/PatCalo versus Pout/Pin (Before Cuts)", 100, 0., 5., 100, 0., 5.);
163 new TH2F(
"MapCor4NoCuts",
"Correlation E25/PatCalo versus E25/highestP (Before Cuts)", 100, 0., 5., 100, 0., 5.);
165 new TH2F(
"MapCor5NoCuts",
"Correlation E25/Pcalo versus Pcalo/Pout (Before Cuts)", 100, 0., 5., 100, 0., 5.);
167 new TH2F(
"MapCor6NoCuts",
"Correlation Pout/Pin versus E25/Pin (Before Cuts)", 100, 0., 5., 100, 0., 5.);
169 new TH2F(
"MapCor7NoCuts",
"Correlation Pout/Pin versus Pcalo/Pout (Before Cuts)", 100, 0., 5., 100, 0., 5.);
171 new TH2F(
"MapCor8NoCuts",
"Correlation E25/Pin versus Pcalo/Pout (Before Cuts)", 100, 0., 5., 100, 0., 5.);
173 new TH2F(
"MapCor9NoCuts",
"Correlation E25/Pcalo versus Eseed/Pout (Before Cuts)", 100, 0., 5., 100, 0., 5.);
175 new TH2F(
"MapCor10NoCuts",
"Correlation Eseed/Pout versus Pout/Pin (Before Cuts)", 100, 0., 5., 100, 0., 5.);
177 new TH2F(
"MapCor11NoCuts",
"Correlation Eseed/Pout versus E25/Pin (Before Cuts)", 100, 0., 5., 100, 0., 5.);
181 Error1NoCuts =
new TH1F(
"Eror1NoCuts",
"DeltaP/Pin (Before Cuts)", 800, -1.0, 1.0);
182 Error2NoCuts =
new TH1F(
"Error2NoCuts",
"DeltaP/Pout (Before Cuts)", 800, -1.0, 1.0);
183 Error3NoCuts =
new TH1F(
"Error3NoCuts",
"DeltaP/Pcalo (Before Cuts)", 800, -1.0, 1.0);
184 eSeedOverPout2NoCuts =
new TH1F(
"eSeedOverPout2NoCuts",
"eSeedOverpout (No Supercluster, Before Cuts)", 600, 0., 4.);
185 hadOverEmNoCuts =
new TH1F(
"hadOverEmNoCuts",
"Had/EM distribution (Before Cuts)", 600, -2., 2.);
189 new TH2F(
"MapCor1ESeed",
"Correlation E25/Pcalo versus E25/Pin (after Eseed/Pout cut)", 100, 0., 5., 100, 0., 5.);
191 new TH2F(
"MapCor2ESeed",
"Correlation E25/Pcalo versus E/P (after Eseed/Pout cut)", 100, 0., 5., 100, 0., 5.);
193 "MapCor3ESeed",
"Correlation E25/Pcalo versus Pout/Pin (after Eseed/Pout cut)", 100, 0., 5., 100, 0., 5.);
195 "MapCor4ESeed",
"Correlation E25/Pcalo versus E25/highestP (after Eseed/Pout cut)", 100, 0., 5., 100, 0., 5.);
197 "MapCor5ESeed",
"Correlation E25/Pcalo versus Pcalo/Pout (after Eseed/Pout cut)", 100, 0., 5., 100, 0., 5.);
199 new TH2F(
"MapCor6ESeed",
"Correlation Pout/Pin versus E25/Pin (after Eseed/Pout cut)", 100, 0., 5., 100, 0., 5.);
201 "MapCor7ESeed",
"Correlation Pout/Pin versus Pcalo/Pout (after Eseed/Pout cut)", 100, 0., 5., 100, 0., 5.);
203 "MapCor8ESeed",
"Correlation E25/Pin versus Pcalo/Pout (after Eseed/Pout cut)", 100, 0., 5., 100, 0., 5.);
205 "MapCor9ESeed",
"Correlation E25/Pcalo versus Eseed/Pout (after Eseed/Pout cut)", 100, 0., 5., 100, 0., 5.);
207 "MapCor10ESeed",
"Correlation Eseed/Pout versus Pout/Pin (after Eseed/Pout cut)", 100, 0., 5., 100, 0., 5.);
209 "MapCor11ESeed",
"Correlation Eseed/Pout versus E25/Pin (after Eseed/Pout cut)", 100, 0., 5., 100, 0., 5.);
212 new TH1F(
"eSeedOverPout2ESeed",
"eSeedOverpout (No Supercluster, after Eseed/Pout cut)", 600, 0., 4.);
214 hadOverEmESeed =
new TH1F(
"hadOverEmESeed",
"Had/EM distribution (after Eseed/Pout cut)", 600, -2., 2.);
217 GeneralMap =
new TH2F(
"GeneralMap",
"Map without any cuts", 173, -86, 86, 362, 0, 361);
218 GeneralMapEndCapMinus =
new TH2F(
"GeneralMapEndCapMinus",
"Map without any cuts", 100, 0, 100, 100, 0, 100);
219 GeneralMapEndCapPlus =
new TH2F(
"GeneralMapEndCapPlus",
"Map without any cuts", 100, 0, 100, 100, 0, 100);
220 GeneralMapBeforePt =
new TH2F(
"GeneralMapBeforePt",
"Map without any cuts", 173, -86, 86, 362, 0, 361);
222 new TH2F(
"GeneralMapEndCapMinusBeforePt",
"Map without any cuts", 100, 0, 100, 100, 0, 100);
224 new TH2F(
"GeneralMapEndCapPlusBeforePt",
"Map without any cuts", 100, 0, 100, 100, 0, 100);
241 <<
" Should be either L3, HH or HHReg. Abort! " << std::endl;
263 time_t cpu_time_used;
267 int nIterations = 10;
283 std::cout <<
" Calibration not run due to problem in Algo Choice..." << std::endl;
304 std::map<EBDetId, float> OldCoeff;
323 std::map<EEDetId, float> OldCoeffEndCap;
336 std::map<DetId, float>::const_iterator itmap;
338 const DetId Id(itmap->first);
340 const EBDetId IChannelDetId(itmap->first);
357 calibs->Fill(itmap->second);
372 write_calibrations.
writeLine(IChannelDetId, itmap->second);
378 const EEDetId IChannelDetId(itmap->first);
386 if (IChannelDetId.
zside() < 0) {
401 write_calibrations.
writeLine(IChannelDetId, itmap->second);
535 std::cout <<
"************* STATISTICS **************" << std::endl;
537 std::cout <<
"Timing info:" << std::endl;
538 std::cout <<
"CPU time usage -- calibrating: " << cpu_time_used <<
" sec." << std::endl;
548 double currEnergy = 0.;
551 for (std::vector<DetId>::const_iterator idsIt = v1.begin(); idsIt != v1.end(); ++idsIt) {
552 if (idsIt->subdetId() == 1) {
554 itrechit = EBhits->
find(*idsIt);
555 if (itrechit == EBhits->
end()) {
556 std::cout <<
"ElectronCalibration::findMaxHit: rechit not found! " << (
EBDetId)(*idsIt) << std::endl;
559 if (itrechit->energy() > currEnergy) {
560 currEnergy = itrechit->energy();
565 itrechit = EEhits->
find(*idsIt);
566 if (itrechit == EEhits->
end()) {
567 std::cout <<
"ElectronCalibration::findMaxHit: rechit not found! idsIt = " << (
EEDetId)(*idsIt) << std::endl;
571 if (itrechit->energy() > currEnergy) {
572 currEnergy = itrechit->energy();
583 if (crystal_ix < 1 || crystal_ix > 100 || crystal_iy < 1 || crystal_iy > 100 ||
abs(iz) != 1) {
586 if ((crystal_ix >= 1 && crystal_ix <= 3 && (crystal_iy <= 40 || crystal_iy > 60)) ||
587 (crystal_ix >= 4 && crystal_ix <= 5 && (crystal_iy <= 35 || crystal_iy > 65)) ||
588 (crystal_ix >= 6 && crystal_ix <= 8 && (crystal_iy <= 25 || crystal_iy > 75)) ||
589 (crystal_ix >= 9 && crystal_ix <= 13 && (crystal_iy <= 20 || crystal_iy > 80)) ||
590 (crystal_ix >= 14 && crystal_ix <= 15 && (crystal_iy <= 15 || crystal_iy > 85)) ||
591 (crystal_ix >= 16 && crystal_ix <= 20 && (crystal_iy <= 13 || crystal_iy > 87)) ||
592 (crystal_ix >= 21 && crystal_ix <= 25 && (crystal_iy <= 8 || crystal_iy > 92)) ||
593 (crystal_ix >= 26 && crystal_ix <= 35 && (crystal_iy <= 5 || crystal_iy > 95)) ||
594 (crystal_ix >= 36 && crystal_ix <= 39 && (crystal_iy <= 3 || crystal_iy > 97)) ||
595 (crystal_ix >= 98 && crystal_ix <= 100 && (crystal_iy <= 40 || crystal_iy > 60)) ||
596 (crystal_ix >= 96 && crystal_ix <= 97 && (crystal_iy <= 35 || crystal_iy > 65)) ||
597 (crystal_ix >= 93 && crystal_ix <= 95 && (crystal_iy <= 25 || crystal_iy > 75)) ||
598 (crystal_ix >= 88 && crystal_ix <= 92 && (crystal_iy <= 20 || crystal_iy > 80)) ||
599 (crystal_ix >= 86 && crystal_ix <= 87 && (crystal_iy <= 15 || crystal_iy > 85)) ||
600 (crystal_ix >= 81 && crystal_ix <= 85 && (crystal_iy <= 13 || crystal_iy > 87)) ||
601 (crystal_ix >= 76 && crystal_ix <= 80 && (crystal_iy <= 8 || crystal_iy > 92)) ||
602 (crystal_ix >= 66 && crystal_ix <= 75 && (crystal_iy <= 5 || crystal_iy > 95)) ||
603 (crystal_ix >= 62 && crystal_ix <= 65 && (crystal_iy <= 3 || crystal_iy > 97)) ||
604 ((crystal_ix == 40 || crystal_ix == 61) &&
605 ((crystal_iy >= 46 && crystal_iy <= 55) || crystal_iy <= 3 || crystal_iy > 97)) ||
606 ((crystal_ix == 41 || crystal_ix == 60) && crystal_iy >= 44 && crystal_iy <= 57) ||
607 ((crystal_ix == 42 || crystal_ix == 59) && crystal_iy >= 43 && crystal_iy <= 58) ||
608 ((crystal_ix == 43 || crystal_ix == 58) && crystal_iy >= 42 && crystal_iy <= 59) ||
609 ((crystal_ix == 44 || crystal_ix == 45 || crystal_ix == 57 || crystal_ix == 56) && crystal_iy >= 41 &&
611 (crystal_ix >= 46 && crystal_ix <= 55 && crystal_iy >= 40 && crystal_iy <= 61)) {
627 std::cerr <<
"Error! can't get the product EBRecHitCollection: " << std::endl;
636 std::cerr <<
"Error! can't get the product EERecHitCollection: " << std::endl;
644 std::cerr <<
"Error! can't get the product ElectronCollection: " << std::endl;
652 if (!EBhits || !EEhits)
674 float highestElePt = 0.;
677 if (fabs(eleIt->eta()) > 2.4)
681 if (eleIt->pt() > highestElePt) {
682 highestElePt = eleIt->pt();
683 highPtElectron = *eleIt;
696 std::vector<DetId> v1;
701 v1.push_back(idsIt->first);
706 maxHitId =
findMaxHit(v1, (EBhits), (EEhits));
710 if (maxHitId.
null()) {
719 maxCC_Eta = ((
EEDetId)maxHitId).ix();
720 maxCC_Phi = ((
EEDetId)maxHitId).iy();
721 Zside = ((
EEDetId)maxHitId).zside();
724 maxCC_Eta = ((
EBDetId)maxHitId).ieta();
725 maxCC_Phi = ((
EBDetId)maxHitId).iphi();
736 std::vector<float>
energy;
737 float energy3x3 = 0.;
738 float energy5x5 = 0.;
745 std::vector<DetId> S9aroundMax = topology->
getWindow(maxHitId, 3, 3);
751 if (S9aroundMax.size() != 9)
758 if (NxNaroundMax[icry].subdetId() ==
EcalBarrel) {
760 itrechit = EBhits->
find(NxNaroundMax[icry]);
761 if (itrechit == EBhits->
end()) {
769 std::cout <<
" nan energy " << std::endl;
772 energy.push_back(itrechit->energy());
773 energy5x5 += itrechit->energy();
776 for (
int tt = 0;
tt < 9;
tt++) {
777 if (NxNaroundMax[icry] == S9aroundMax[
tt])
778 energy3x3 += itrechit->energy();
783 itrechit = EEhits->
find(NxNaroundMax[icry]);
785 if (itrechit == EEhits->
end()) {
794 std::cout <<
" nan energy " << std::endl;
797 energy.push_back(itrechit->energy());
798 energy5x5 += itrechit->energy();
801 for (
int tt = 0;
tt < 9;
tt++) {
802 if (NxNaroundMax[icry] == S9aroundMax[
tt])
803 energy3x3 += itrechit->energy();
825 if (highestElePt <
ElePt_)
846 float UncorrectedPatCalo =
873 MapCor1NoCuts->Fill(energy5x5 / UncorrectedPatCalo, energy5x5 / Ptrack_in);
875 MapCor3NoCuts->Fill(energy5x5 / UncorrectedPatCalo, Ptrack_out / Ptrack_in);
876 MapCor4NoCuts->Fill(energy5x5 / UncorrectedPatCalo, energy5x5 / highPtElectron.
p());
877 MapCor5NoCuts->Fill(energy5x5 / UncorrectedPatCalo, UncorrectedPatCalo / Ptrack_out);
878 MapCor6NoCuts->Fill(Ptrack_out / Ptrack_in, energy5x5 / Ptrack_in);
879 MapCor7NoCuts->Fill(Ptrack_out / Ptrack_in, UncorrectedPatCalo / Ptrack_out);
880 MapCor8NoCuts->Fill(energy5x5 / Ptrack_in, UncorrectedPatCalo / Ptrack_out);
895 if ((energy3x3 / energy5x5) <
cut1_)
897 if ((Ptrack_out / Ptrack_in) <
cut2_ || (Ptrack_out / Ptrack_in) >
cut3_)
972 E25oPvsEta->Fill(maxCC_Eta, energy5x5 / UncorrectedPatCalo);
974 float Radius =
sqrt((maxCC_Eta) * (maxCC_Eta) + (maxCC_Phi) * (maxCC_Phi));
982 e25->Fill(energy5x5);
985 trP->Fill(UncorrectedPatCalo);
990 E25oP->Fill(energy5x5 / UncorrectedPatCalo);
993 Map->Fill(maxCC_Eta, maxCC_Phi);
1009 MapCor1->Fill(energy5x5 / UncorrectedPatCalo, energy5x5 / Ptrack_in);
1011 MapCor3->Fill(energy5x5 / UncorrectedPatCalo, Ptrack_out / Ptrack_in);
1012 MapCor4->Fill(energy5x5 / UncorrectedPatCalo, energy5x5 / highPtElectron.
p());
1013 MapCor5->Fill(energy5x5 / UncorrectedPatCalo, UncorrectedPatCalo / Ptrack_out);
1014 MapCor6->Fill(Ptrack_out / Ptrack_in, energy5x5 / Ptrack_in);
1015 MapCor7->Fill(Ptrack_out / Ptrack_in, UncorrectedPatCalo / Ptrack_out);
1016 MapCor8->Fill(energy5x5 / Ptrack_in, UncorrectedPatCalo / Ptrack_out);
1021 PinMinPout->Fill((Ptrack_in - Ptrack_out) / Ptrack_in);