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