00001 #include "Calibration/IsolatedParticles/plugins/ElectronStudy.h"
00002 #include "Calibration/IsolatedParticles/interface/CaloPropagateTrack.h"
00003 #include "Calibration/IsolatedParticles/interface/eECALMatrix.h"
00004 #include "FWCore/Utilities/interface/Exception.h"
00005 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
00006 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
00007 #include "Geometry/CaloEventSetup/interface/CaloTopologyRecord.h"
00008 #include "Geometry/CaloTopology/interface/CaloSubdetectorTopology.h"
00009 #include "Geometry/CaloTopology/interface/CaloTopology.h"
00010 #include "MagneticField/Engine/interface/MagneticField.h"
00011 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
00012
00013 #include <algorithm>
00014
00015 ElectronStudy::ElectronStudy(const edm::ParameterSet& ps) {
00016
00017 sourceLabel = ps.getUntrackedParameter<std::string>("SourceLabel","generator");
00018 g4Label = ps.getUntrackedParameter<std::string>("ModuleLabel","g4SimHits");
00019 hitLabEB= ps.getUntrackedParameter<std::string>("EBCollection","EcalHitsEB");
00020 hitLabEE= ps.getUntrackedParameter<std::string>("EECollection","EcalHitsEE");
00021 verbose = ps.getUntrackedParameter<int>("Verbosity",0);
00022 edm::LogInfo("ElectronStudy") << "Module Label: " << g4Label << " Hits: "
00023 << hitLabEB << ", " << hitLabEE;
00024
00025 double tempPBins[NPBins+1] = { 0.0, 10.0, 20.0, 40.0, 60.0,
00026 100.0, 500.0, 1000.0, 10000.0};
00027 double tempEta[NEtaBins+1] = {0.0, 1.2, 1.6, 3.0};
00028
00029 for (int i=0; i<NPBins+1; i++) pBins[i] = tempPBins[i];
00030 for(int i=0; i<NEtaBins+1; i++) etaBins[i] = tempEta[i];
00031
00032 edm::Service<TFileService> tfile;
00033 if ( !tfile.isAvailable() ) {
00034 edm::LogInfo("ElectronStudy") << "TFileService unavailable: no histograms";
00035 histos = false;
00036 } else {
00037 char name[20], title[200], cpbin[30], cebin[30];
00038 histos = true;
00039 for (unsigned int i=0; i<NPBins+1; ++i) {
00040 if (i == 0) sprintf (cpbin, " All p");
00041 else sprintf (cpbin, " p (%6.0f:%6.0f)", pBins[i-1], pBins[i]);
00042 for (unsigned int j=0; j<NEtaBins+1; ++j) {
00043 if (j == 0) sprintf (cebin, " All #eta");
00044 else sprintf (cebin, " #eta (%4.1f:%4.1f)", etaBins[j-1], etaBins[j]);
00045 sprintf (name, "R1%d%d", i, j);
00046 sprintf (title,"E1/E9 for %s%s", cpbin, cebin);
00047 histoR1[i][j] = tfile->make<TH1F>(name, title, 100, 0., 2.);
00048 histoR1[i][j]->GetXaxis()->SetTitle(title);
00049 histoR1[i][j]->GetYaxis()->SetTitle("Tracks");
00050 sprintf (name, "R2%d%d", i, j);
00051 sprintf (title,"E1/E25 for %s%s", cpbin, cebin);
00052 histoR2[i][j] = tfile->make<TH1F>(name, title, 100, 0., 2.);
00053 histoR2[i][j]->GetXaxis()->SetTitle(title);
00054 histoR2[i][j]->GetYaxis()->SetTitle("Tracks");
00055 sprintf (name, "R3%d%d", i, j);
00056 sprintf (title,"E9/E25 for %s%s", cpbin, cebin);
00057 histoR3[i][j] = tfile->make<TH1F>(name, title, 100, 0., 2.);
00058 histoR3[i][j]->GetXaxis()->SetTitle(title);
00059 histoR3[i][j]->GetYaxis()->SetTitle("Tracks");
00060 sprintf (name, "E1x1%d%d", i, j);
00061 sprintf (title,"E1/P for %s%s", cpbin, cebin);
00062 histoE1x1[i][j] = tfile->make<TH1F>(name, title, 100, 0., 2.);
00063 histoE1x1[i][j]->GetXaxis()->SetTitle(title);
00064 histoE1x1[i][j]->GetYaxis()->SetTitle("Tracks");
00065 sprintf (name, "E3x3%d%d", i, j);
00066 sprintf (title,"E9/P for %s%s", cpbin, cebin);
00067 histoE3x3[i][j] = tfile->make<TH1F>(name, title, 100, 0., 2.);
00068 histoE3x3[i][j]->GetXaxis()->SetTitle(title);
00069 histoE3x3[i][j]->GetYaxis()->SetTitle("Tracks");
00070 sprintf (name, "E5x5%d%d", i, j);
00071 sprintf (title,"E25/P for %s%s", cpbin, cebin);
00072 histoE5x5[i][j] = tfile->make<TH1F>(name, title, 100, 0., 2.);
00073 histoE5x5[i][j]->GetXaxis()->SetTitle(title);
00074 histoE5x5[i][j]->GetYaxis()->SetTitle("Tracks");
00075 }
00076 }
00077 }
00078 }
00079
00080 void ElectronStudy::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
00081
00082 if (verbose > 1) std::cout << "Run = " << iEvent.id().run() << " Event = "
00083 << iEvent.id().event() << std::endl;
00084
00085
00086 edm::ESHandle<MagneticField> bFieldH;
00087 iSetup.get<IdealMagneticFieldRecord>().get(bFieldH);
00088 const MagneticField* bField = bFieldH.product();
00089
00090 edm::ESHandle<CaloGeometry> pG;
00091 iSetup.get<CaloGeometryRecord>().get(pG);
00092 const CaloGeometry* geo = pG.product();
00093
00094 edm::ESHandle<CaloTopology> theCaloTopology;
00095 iSetup.get<CaloTopologyRecord>().get(theCaloTopology);
00096 const CaloTopology* caloTopology = theCaloTopology.product();
00097
00098
00099 edm::Handle<edm::PCaloHitContainer> caloHitEB;
00100 iEvent.getByLabel(g4Label,hitLabEB,caloHitEB);
00101
00102
00103 edm::Handle<edm::PCaloHitContainer> caloHitEE;
00104 iEvent.getByLabel(g4Label,hitLabEE,caloHitEE);
00105
00106
00107 edm::Handle<edm::SimTrackContainer> SimTk;
00108 iEvent.getByLabel(g4Label, SimTk);
00109
00110
00111 edm::Handle<edm::SimVertexContainer> SimVtx;
00112 iEvent.getByLabel(g4Label, SimVtx);
00113
00114 if (verbose>0)
00115 std::cout << "ElectronStudy: hits valid[EB]: " << caloHitEB.isValid()
00116 << " valid[EE]: " << caloHitEE.isValid() << std::endl;
00117
00118 if (caloHitEB.isValid() && caloHitEE.isValid()) {
00119 unsigned int indx;
00120 if (verbose>2) {
00121 edm::PCaloHitContainer::const_iterator ihit;
00122 for (ihit=caloHitEB->begin(),indx=0; ihit!=caloHitEB->end(); ihit++,indx++) {
00123 EBDetId id = ihit->id();
00124 std::cout << "Hit[" << indx << "] " << id << " E " << ihit->energy()
00125 << " T " << ihit->time() << std::endl;
00126 }
00127 for (ihit=caloHitEE->begin(),indx=0; ihit!=caloHitEE->end(); ihit++,indx++) {
00128 EEDetId id = ihit->id();
00129 std::cout << "Hit[" << indx << "] " << id << " E " << ihit->energy()
00130 << " T " << ihit->time() << std::endl;
00131 }
00132 }
00133 edm::SimTrackContainer::const_iterator simTrkItr=SimTk->begin();
00134 for (indx=0; simTrkItr!= SimTk->end(); simTrkItr++,indx++) {
00135 if (verbose>0) std::cout << "ElectronStudy: Track[" << indx << "] ID "
00136 << simTrkItr->trackId() << " type "
00137 << simTrkItr->type() << " charge "
00138 << simTrkItr->charge() << " p "
00139 << simTrkItr->momentum()<< " Generator Index "
00140 << simTrkItr->genpartIndex() << " vertex "
00141 << simTrkItr->vertIndex() << std::endl;
00142 if (std::abs(simTrkItr->type()) == 11 && simTrkItr->vertIndex() != -1) {
00143 int thisTrk = simTrkItr->trackId();
00144 spr::propagatedTrackDirection trkD = spr::propagateCALO(thisTrk, SimTk, SimVtx, geo, bField, (verbose>1));
00145 if (trkD.okECAL) {
00146 const DetId isoCell = trkD.detIdECAL;
00147 double e1x1 = spr::eECALmatrix(isoCell, caloHitEB, caloHitEE, geo, caloTopology, 0, 0, -100.0, -100.0,-500.0, 500.0, (verbose>2));
00148 double e3x3 = spr::eECALmatrix(isoCell, caloHitEB, caloHitEE, geo, caloTopology, 1, 1, -100.0, -100.0,-500.0, 500.0, (verbose>2));
00149 double e5x5 = spr::eECALmatrix(isoCell, caloHitEB, caloHitEE, geo, caloTopology, 2, 2, -100.0, -100.0,-500.0, 500.0, (verbose>2));
00150 double p = simTrkItr->momentum().P();
00151 double eta = std::abs(simTrkItr->momentum().eta());
00152 int etaBin=-1, momBin=-1;
00153 for (int ieta=0; ieta<NEtaBins; ieta++) {
00154 if (eta>etaBins[ieta] && eta<etaBins[ieta+1] ) etaBin = ieta+1;
00155 }
00156 for (int ipt=0; ipt<NPBins; ipt++) {
00157 if (p>pBins[ipt] && p<pBins[ipt+1] ) momBin = ipt+1;
00158 }
00159 double r1=-1, r2=-1, r3=-1;
00160 if (e3x3 > 0) r1 = e1x1/e3x3;
00161 if (e5x5 > 0) {r2 = e1x1/e5x5; r3 = e3x3/e5x5;}
00162 if (verbose>0) {
00163 std::cout << "ElectronStudy: p " << p << " [" << momBin << "] eta "
00164 << eta << " [" << etaBin << "]";
00165 if (isoCell.subdetId() == EcalBarrel) {
00166 EBDetId id = isoCell;
00167 std::cout << " Cell 0x" << std::hex << isoCell() << std::dec
00168 << " " << id;
00169 } else if (isoCell.subdetId() == EcalEndcap) {
00170 EEDetId id = isoCell;
00171 std::cout << " Cell 0x" << std::hex << isoCell() << std::dec
00172 << " " << id;
00173 } else {
00174 std::cout << " Cell 0x" << std::hex << isoCell() << std::dec
00175 << " Unknown Type";
00176 }
00177 std::cout << " e1x1 " << e1x1 << "|" << r1 << "|" << r2 << " e3x3 "
00178 << e3x3 << "|" << r3 << " e5x5 " << e5x5 << std::endl;
00179 }
00180 if (histos) {
00181 histoR1[0][0]->Fill(r1);
00182 histoR2[0][0]->Fill(r2);
00183 histoR3[0][0]->Fill(r3);
00184 histoE1x1[0][0]->Fill(e1x1/p);
00185 histoE3x3[0][0]->Fill(e3x3/p);
00186 histoE5x5[0][0]->Fill(e5x5/p);
00187 if (momBin>0) {
00188 histoR1[momBin][0]->Fill(r1);
00189 histoR2[momBin][0]->Fill(r2);
00190 histoR3[momBin][0]->Fill(r3);
00191 histoE1x1[momBin][0]->Fill(e1x1/p);
00192 histoE3x3[momBin][0]->Fill(e3x3/p);
00193 histoE5x5[momBin][0]->Fill(e5x5/p);
00194 }
00195 if (etaBin>0) {
00196 histoR1[0][etaBin]->Fill(r1);
00197 histoR2[0][etaBin]->Fill(r2);
00198 histoR3[0][etaBin]->Fill(r3);
00199 histoE1x1[0][etaBin]->Fill(e1x1/p);
00200 histoE3x3[0][etaBin]->Fill(e3x3/p);
00201 histoE5x5[0][etaBin]->Fill(e5x5/p);
00202 if (momBin>0) {
00203 histoR1[momBin][etaBin]->Fill(r1);
00204 histoR2[momBin][etaBin]->Fill(r2);
00205 histoR3[momBin][etaBin]->Fill(r3);
00206 histoE1x1[momBin][etaBin]->Fill(e1x1/p);
00207 histoE3x3[momBin][etaBin]->Fill(e3x3/p);
00208 histoE5x5[momBin][etaBin]->Fill(e5x5/p);
00209 }
00210 }
00211 }
00212 }
00213 }
00214 }
00215 }
00216
00217 }
00218
00219
00220 DEFINE_FWK_MODULE(ElectronStudy);
1.7.3