00001
00002
00003
00004
00005
00006 #include "MCCheckOut/CosmicAna.h"
00007 #include <vector>
00008 #include <cmath>
00009 #include "TH1F.h"
00010 #include "TH2F.h"
00011 #include "Config/Config.h"
00012 #include "EventDataModel/edm.h"
00013 #include "Geometry/geo.h"
00014 #include "JobControl/jobc.h"
00015 #include "Simulation/PartToHitTable.h"
00016 #include "Simulation/ParticleList.h"
00017 #include "Simulation/FLSHitList.h"
00018 using namespace mcchk;
00019 MODULE_DECL(mcchk::CosmicAna);
00020
00021
00022
00023 CosmicAna::CosmicAna(const char* version) :
00024 jobc::Module("mcchk::CosmicAna")
00025 {
00026 this->SetCfgVersion(version);
00027
00028 fNhitHisto = new TH1F("fNhitHisto",
00029 ";hits/20 usec;spills",200,0.0,20000.0);
00030
00031 fDminHisto0 = new TH1F("fDminHisto0",";d (cm);",100,0.0,500.E2);
00032 fDminHisto1 = new TH1F("fDminHisto1",";d (cm);",100,0.0,500.E2);
00033
00034 fPhotonAngles = new TH2F("fPhotonAngles",";x;y",
00035 36,-180.0,180.0,50,-1.0,1.0);
00036 fPhotonAnglesLo = new TH2F("fPhotonAnglesLo",";x;y",
00037 36,-180.0,180.0,50,-1.0,1.0);
00038 fPhotonAnglesMi = new TH2F("fPhotonAnglesMi",";x;y",
00039 36,-180.0,180.0,50,-1.0,1.0);
00040 fPhotonAnglesHi = new TH2F("fPhotonAnglesHi",";x;y",
00041 36,-180.0,180.0,50,-1.0,1.0);
00042 fPhotonCosQ = new TH1F("fPhotonCosQ",";cosQ;tracks",50,-1.0,1.0);
00043 fPhotonEnergy = new TH1F("fPhotonEnergy",";E (GeV)",5000,0.0,1000.0);
00044
00045 fElectronAngles = new TH2F("fElectronAngles",";x;y",
00046 36,-180.0,180.0,50,-1.0,1.0);
00047 fElectronAnglesLo = new TH2F("fElectronAnglesLo",";x;y",
00048 36,-180.0,180.0,50,-1.0,1.0);
00049 fElectronAnglesMi = new TH2F("fElectronAnglesMi",";x;y",
00050 36,-180.0,180.0,50,-1.0,1.0);
00051 fElectronAnglesHi = new TH2F("fElectronAnglesHi",";x;y",
00052 36,-180.0,180.0,50,-1.0,1.0);
00053 fElectronCosQ = new TH1F("fElectronCosQ",";cosQ;tracks",50,-1.0,1.0);
00054 fElectronEnergy = new TH1F("fElectronEnergy",";E (GeV)",5000,0.0,1000.0);
00055
00056 fMuonAngles = new TH2F("fMuonAngles",";x;y",
00057 36,-180.0,180.0,50,-1.0,1.0);
00058 fMuonAnglesLo = new TH2F("fMuonAnglesLo",";x;y",
00059 36,-180.0,180.0,50,-1.0,1.0);
00060 fMuonAnglesMi = new TH2F("fMuonAnglesMi",";x;y",
00061 36,-180.0,180.0,50,-1.0,1.0);
00062 fMuonAnglesHi = new TH2F("fMuonAnglesHi",";x;y",
00063 36,-180.0,180.0,50,-1.0,1.0);
00064 fMuonCosQ = new TH1F("fMuonCosQ",";cosQ;tracks",50,-1.0,1.0);
00065 fMuonEnergy = new TH1F("fMuonEnergy",";E (GeV)",5000,0.0,1000.0);
00066 }
00067
00068
00069
00070 CosmicAna::~CosmicAna() { }
00071
00072
00073
00074 void CosmicAna::Update(const cfg::Config& ) { }
00075
00076
00077
00078 jobc::Result CosmicAna::Ana(const edm::EventHandle& evt)
00079 {
00080 std::vector<const sim::ParticleList*> partlist;
00081 try { evt.DetSim().Get(".",partlist); }
00082 catch (...) { std::cerr << "No particle list!" << std::endl; }
00083
00084 std::vector<const sim::FLSHitList*> hitlist;
00085 try { evt.DetSim().Get(".",hitlist); }
00086 catch (...) { std::cerr << "No FLS hit list!" << std::endl; }
00087
00088 geo::Geometry& geom = geo::Geometry::Instance(evt.Header().Run(),1);
00089
00090 const std::vector<sim::Particle>& particle = partlist[0]->fParticles;
00091 const std::vector<sim::FLSHit>& hits = hitlist[0]->fHits;
00092 fNhitHisto->Fill((float)hits.size());
00093
00094
00095 sim::PartToHitTable p2h(*partlist[0], *hitlist[0]);
00096
00097 double electronE = -1.0;
00098 double photonE = -1.0;
00099 for (unsigned int i=0; i<particle.size(); ++i) {
00100 const TLorentzVector& v4 = particle[i].fVtx;
00101 const TLorentzVector& p4 = particle[i].fP;
00102 double x0[3] = {v4.X(), v4.Y(), v4.Z()};
00103 double dx[3] = {p4.Px(), p4.Py(), p4.Pz()};
00104 double x1[3] = {0.0, 0.0, 0.5*geom.DetLength()};
00105 double x2[3];
00106 double d = geo::ClosestApproach(x1, x0, dx, x2);
00107 if (p2h[i].size()==0) {fDminHisto0->Fill(d); }
00108 else {fDminHisto1->Fill(d); }
00109 if (p2h[i].size()>0 && d>100.E2) {
00110 std::cout <<
00111 "dmin.gt.100 " <<
00112 evt.Header().Run() << " " << evt.Header().Event() << " " <<
00113 particle[i].fPdgCode << " " <<
00114 v4.X() << " " << v4.Y() << " " << v4.Z() << " " <<
00115 p4.Px() << " " << p4.Py() << " " << p4.Pz() << " " <<
00116 d << " " << p2h[i].size() << std::endl;
00117 }
00118
00119 TH1F* hCosQ = 0;
00120 TH2F* hAngles = 0;
00121 TH2F* hAnglesLo = 0;
00122 TH2F* hAnglesMi = 0;
00123 TH2F* hAnglesHi = 0;
00124 TH1F* hEnergy = 0;
00125 if (abs(particle[i].fPdgCode)==13) {
00126 hCosQ = fMuonCosQ;
00127 hAngles = fMuonAngles;
00128 hAnglesLo = fMuonAnglesLo;
00129 hAnglesMi = fMuonAnglesMi;
00130 hAnglesHi = fMuonAnglesHi;
00131 hEnergy = fMuonEnergy;
00132 }
00133 else if (abs(particle[i].fPdgCode)==22) {
00134 hCosQ = fPhotonCosQ;
00135 hAngles = fPhotonAngles;
00136 hAnglesLo = fPhotonAnglesLo;
00137 hAnglesMi = fPhotonAnglesMi;
00138 hAnglesHi = fPhotonAnglesHi;
00139 hEnergy = fPhotonEnergy;
00140 }
00141 else if (abs(particle[i].fPdgCode)==11) {
00142 hCosQ = fElectronCosQ;
00143 hAngles = fElectronAngles;
00144 hAnglesLo = fElectronAnglesLo;
00145 hAnglesMi = fElectronAnglesMi;
00146 hAnglesHi = fElectronAnglesHi;
00147 hEnergy = fElectronEnergy;
00148 }
00149 if (hCosQ!=0) {
00150 if (p2h[i].size()>0) {
00151 double cosq = -p4.Py()/p4.P();
00152 double phi = atan2(p4.Pz(),p4.Px());
00153 phi *= 180/M_PI;
00154 hCosQ->Fill(cosq);
00155 hAngles->Fill(phi,cosq);
00156 if (p4.E()<1.0) hAnglesLo->Fill(phi,cosq);
00157 else if (p4.E()<10.0) hAnglesMi->Fill(phi,cosq);
00158 else hAnglesHi->Fill(phi,cosq);
00159 hEnergy->Fill(p4.E());
00160
00161
00162 if (p4.E()>0.500 && p4.E()<5.0) {
00163 if (abs(particle[i].fPdgCode)==11) electronE = p4.E();
00164 if (abs(particle[i].fPdgCode)==22) photonE = p4.E();
00165 }
00166 }
00167 }
00168 }
00169 if (electronE>0.0 || photonE>0.0) {
00170 std::cout << "Run/Event " << evt.Header().Run() << " "
00171 << evt.Header().Event() << " "
00172 << electronE << " " << photonE << std::endl;
00173 }
00174 return jobc::kPassed;
00175 }
00176
