NOvA Offline: mcchk::CosmicAna Class Reference

00095 { 
00096   std::vector<const rawdata::RawDigit*> digits;
00097   try         { evt.Raw().Get(".",digits); }
00098   catch (...) { 
00099     if(norawcnt < 5){
00100       ++norawcnt;
00101       std::cerr << "No raw digit list!" << std::endl; 
00102     }
00103     return jobc::kFailed;
00104   }
00105     
00106   std::vector<const sim::ParticleList*> partlist;
00107   try         { evt.DetSim().Get(".",partlist); }
00108   catch (...) { 
00109     if(nopartcnt < 5){
00110       ++nopartcnt;
00111       std::cerr << "No particle list!" << std::endl; 
00112     }
00113     return jobc::kFailed;
00114   }
00115 
00116   std::vector<const sim::FLSHitList*> hitlist;
00117   try         { evt.DetSim().Get(".",hitlist); }
00118   catch (...) { 
00119     if(noflscnt < 5){
00120       ++noflscnt;
00121       std::cerr << "No FLS hit list!" << std::endl; 
00122     }
00123     return jobc::kFailed;
00124   }
00125     
00126   geo::Geometry& geom = util::EDMUtils::GetGeometry(evt);
00127   cmap::CMap& cmap = util::EDMUtils::GetCMap(evt);
00128 
00129   const std::vector<sim::Particle>& particle = partlist[0]->fParticles;
00130   const std::vector<sim::FLSHit>&   hits     = hitlist[0]->fHits;
00131   fNhitHisto->Fill((float)hits.size());
00132 
00134   for(unsigned int i = 0; i < digits.size(); ++i){
00135     double xyz[3] = {0.};
00136 
00137     // Need to ask channel map for plane and cell number.
00138     geom.Plane(cmap.GetPlane(digits[i])).Cell(cmap.GetCell(digits[i])).GetCenter(xyz);
00139 
00140 //     std::cout << i 
00141 //            << " " << xyz[0]
00142 //            << " " << xyz[1]
00143 //            << " " << xyz[2]
00144 //            << " " << digits[i]->GetPlane()
00145 //            << " " << digits[i]->GetCell()
00146 //            << " " << geom.Plane(digits[i]->GetPlane()).View()
00147 //            << std::endl;
00148   }
00149 
00151   double wpos[3] = {0.};
00152   for(unsigned int i = 0; i < hits.size(); ++i){
00154     double lpos[3] = {0.5*(hits[i].fEntryX+hits[i].fExitX),
00155                       0.5*(hits[i].fEntryY+hits[i].fExitY),
00156                       0.5*(hits[i].fEntryZ+hits[i].fExitZ)};
00157     geom.Plane(hits[i].fPlaneId).Cell(hits[i].fCellId).LocalToWorld(lpos,wpos);
00158 
00159 //     std::cout << evt.Header().Event() << " " << i 
00160 //            << " " << hits[i].fPDG 
00161 //            << " " << hits[i].fTrackId
00162 //            << " " << hits[i].fPlaneId
00163 //            << " " << hits[i].fCellId 
00164 //            << " " << wpos[0] 
00165 //            << " " << wpos[1] 
00166 //            << " " << wpos[2] 
00167 //            << std::endl;
00168 
00169     if(geom.Plane(hits[i].fPlaneId).View() == geo::kY)
00170       fHitsYZ->Fill(wpos[2], wpos[1]);
00171     //fHitsYZ->Fill(hits[i].fPlaneId, hits[i].fCellId);
00172     if(geom.Plane(hits[i].fPlaneId).View() == geo::kX)
00173       fHitsXZ->Fill(wpos[2], wpos[0]);
00174     //fHitsXZ->Fill(hits[i].fPlaneId, hits[i].fCellId);
00175   }
00176 
00177   // Build look up table from a particle to the hits it created
00178   sim::PartToHitTable p2h(*partlist[0], *hitlist[0]);
00179   
00180   double electronE = -1.0;
00181   double photonE   = -1.0;
00182   for (unsigned int i=0; i<particle.size(); ++i) {
00183     const TLorentzVector& v4 = particle[i].fVtx;
00184     const TLorentzVector& p4 = particle[i].fP;
00185     double x0[3] = {v4.X(),  v4.Y(),  v4.Z()};
00186     double dx[3] = {p4.Px(), p4.Py(), p4.Pz()};
00187     double x1[3] = {0.0,     0.0,     0.5*geom.DetLength()};
00188     double x2[3];
00189     double d = geo::ClosestApproach(x1, x0, dx, x2);
00190     if (p2h[i].size()==0) {fDminHisto0->Fill(d); }
00191     else                  {fDminHisto1->Fill(d); }
00192     if (p2h[i].size()>0 && d>100.E2) {
00193       std::cout <<
00194         "dmin.gt.100 " << 
00195         evt.Header().Run() << " " << evt.Header().Event() << " " <<
00196         particle[i].fPdgCode << " " <<
00197         v4.X()  << " " << v4.Y()  << " " << v4.Z()  << " " <<
00198         p4.Px() << " " << p4.Py() << " " << p4.Pz() << " " <<
00199         d << " " << p2h[i].size() << std::endl;
00200     }
00201 
00202     TH1F* hCosQ     = 0;
00203     TH2F* hAngles   = 0;
00204     TH2F* hAnglesLo = 0;
00205     TH2F* hAnglesMi = 0;
00206     TH2F* hAnglesHi = 0;
00207     TH1F* hEnergy   = 0;
00208     if (abs(particle[i].fPdgCode)==13) {
00209       hCosQ     = fMuonCosQ;
00210       hAngles   = fMuonAngles;
00211       hAnglesLo = fMuonAnglesLo;
00212       hAnglesMi = fMuonAnglesMi;
00213       hAnglesHi = fMuonAnglesHi;
00214       hEnergy   = fMuonEnergy;
00215     }
00216     else if (abs(particle[i].fPdgCode)==22) {
00217       hCosQ     = fPhotonCosQ;
00218       hAngles   = fPhotonAngles;
00219       hAnglesLo = fPhotonAnglesLo;
00220       hAnglesMi = fPhotonAnglesMi;
00221       hAnglesHi = fPhotonAnglesHi;
00222       hEnergy   = fPhotonEnergy;
00223     }
00224     else if (abs(particle[i].fPdgCode)==11) {
00225       hCosQ     = fElectronCosQ;
00226       hAngles   = fElectronAngles;
00227       hAnglesLo = fElectronAnglesLo;
00228       hAnglesMi = fElectronAnglesMi;
00229       hAnglesHi = fElectronAnglesHi;
00230       hEnergy   = fElectronEnergy;
00231     }
00232     if (hCosQ!=0) {
00233       if (p2h[i].size()>0) {
00234         double cosq = -p4.Py()/p4.P();
00235         double phi  = atan2(p4.Pz(),p4.Px());
00236         phi *= 180/M_PI;
00237         hCosQ->Fill(cosq);
00238         hAngles->Fill(phi,cosq);
00239         if      (p4.E()<1.0)  hAnglesLo->Fill(phi,cosq);
00240         else if (p4.E()<10.0) hAnglesMi->Fill(phi,cosq);
00241         else                  hAnglesHi->Fill(phi,cosq);
00242         hEnergy->Fill(p4.E());
00243 
00244         fCRInfo.run = evt.Header().Run();
00245         fCRInfo.event = evt.Header().Event();
00246         fCRInfo.pdg = particle[i].fPdgCode;
00247         fCRInfo.dcosX = p4.Px()/p4.P();
00248         fCRInfo.dcosY = p4.Py()/p4.P();
00249         fCRInfo.dcosZ = p4.Pz()/p4.P();
00250         fCRInfo.energy = p4.E();
00251 
00252         std::vector< std::vector<double> > hitsP(geom.GetPlanesByView().size());
00253         unsigned int numCells = geom.Plane(0).Ncells();
00254         if(geom.Plane(1).Ncells() > numCells) numCells = geom.Plane(1).Ncells();
00255         for(unsigned int p = 0; p < hitsP.size(); ++p)
00256           hitsP[p].resize(numCells,0.);
00257 
00258         std::vector<geo::View_t> viewP(geom.GetPlanesByView().size());
00259         
00260         double entryPos[3] = {0.};
00261         double exitPos[3] = {0.};
00262         for(unsigned int h = 0; h < p2h[i].size(); ++h){
00263           if(h == 0){
00264             double lpos[3] = {0.5*(p2h[i][h]->fEntryX+p2h[i][h]->fExitX),
00265                               0.5*(p2h[i][h]->fEntryY+p2h[i][h]->fExitY),
00266                               0.5*(p2h[i][h]->fEntryZ+p2h[i][h]->fExitZ)};
00267             geom.Plane(p2h[i][h]->fPlaneId).Cell(p2h[i][h]->fCellId).LocalToWorld(lpos,entryPos);
00268           }
00269           double lpos[3] = {0.5*(p2h[i][h]->fEntryX+p2h[i][h]->fExitX),
00270                             0.5*(p2h[i][h]->fEntryY+p2h[i][h]->fExitY),
00271                             0.5*(p2h[i][h]->fEntryZ+p2h[i][h]->fExitZ)};
00272           geom.Plane(p2h[i][h]->fPlaneId).Cell(p2h[i][h]->fCellId).LocalToWorld(lpos,exitPos);
00273           
00274           hitsP[p2h[i][h]->fPlaneId][p2h[i][h]->fCellId] += 1.;
00275           viewP[p2h[i][h]->fPlaneId] = geom.Plane(p2h[i][h]->fPlaneId).View();
00276 
00277 //        std::cout << evt.Header().Event() << " " << i << " " << h << " " 
00278 //                  << fCRInfo.pdg << " " << p2h[i][h]->fPlaneId 
00279 //                  << " " << p2h[i][h]->fCellId 
00280 //                  << " " << entryPos[0] 
00281 //                  << " " << entryPos[1] 
00282 //                  << " " << entryPos[2] 
00283 //                  << " " << exitPos[0] 
00284 //                  << " " << exitPos[1] 
00285 //                  << " " << exitPos[2] 
00286 //                  << std::endl;
00287 
00288         }
00289 
00290         double planesX = 0.;
00291         double cellsX = 0.;
00292         double planesY = 0.;
00293         double cellsY = 0.;
00294         for(unsigned int p = 0; p < hitsP.size(); ++p){
00295           double cells = 0.;
00296           for(unsigned int c = 0; c < hitsP[p].size(); ++c){
00297             if(hitsP[p][c] > 0.) cells += 1.;
00298           }
00299            
00300           if(viewP[p] == geo::kX && cells > 0.){
00301             planesX += 1.;
00302             cellsX += cells;
00303           }
00304           if(viewP[p] == geo::kY && cells > 0.){
00305             planesY += 1.;
00306             cellsY += cells;
00307           }
00308 
00309 //        std::cout << "plane " << p << " has " << cells
00310 //                  << " hit cells" << std::endl;
00311 
00312         }//end loop over hit cells
00313 
00314 
00315         if(planesX > 0.) fCRInfo.avCellsX = cellsX/planesX;
00316         else             fCRInfo.avCellsX = 0.;
00317   
00318         if(planesY > 0.) fCRInfo.avCellsY = cellsY/planesY;
00319         else             fCRInfo.avCellsY = 0.;
00320 
00321         fCRInfo.vtxX = entryPos[0];
00322         fCRInfo.vtxY = entryPos[1];
00323         fCRInfo.vtxZ = entryPos[2];
00324         fCRInfo.endX = exitPos[0];
00325         fCRInfo.endY = exitPos[1];
00326         fCRInfo.endZ = exitPos[2];
00327 
00329         if(exitPos[1] > entryPos[1]){
00330           fCRInfo.endX = entryPos[0];
00331           fCRInfo.endY = entryPos[1];
00332           fCRInfo.endZ = entryPos[2];
00333           fCRInfo.vtxX = exitPos[0];
00334           fCRInfo.vtxY = exitPos[1];
00335           fCRInfo.vtxZ = exitPos[2];
00336         }
00337 
00338         //check this calculation
00339         fCRInfo.length = TMath::Sqrt(TMath::Power(exitPos[0]-entryPos[0], 2.)
00340                                      + TMath::Power(exitPos[1]-entryPos[1], 2.)
00341                                      + TMath::Power(exitPos[2]-entryPos[2], 2.));
00342 
00343 
00344         fTree->Fill();
00345 
00346         // Check for interesting particles
00347         if (p4.E()>0.500 && p4.E()<5.0) {
00348           if (abs(particle[i].fPdgCode)==11) electronE = p4.E();
00349           if (abs(particle[i].fPdgCode)==22) photonE   = p4.E();
00350         }
00351       }
00352     } // loop on particles
00353   }
00354   if (electronE>0.0 || photonE>0.0) {
00355     std::cout << "Run/Event " << evt.Header().Run() << " "
00356               << evt.Header().Event() << " " 
00357               << electronE << " " << photonE << std::endl;
00358   }
00359   return jobc::kPassed;
00360 }


Public Member Functions
	CosmicAna (const char *version)
	~CosmicAna ()
void	Update (const cfg::Config &c)
jobc::Result	Ana (const edm::EventHandle &evt)
Private Attributes
TH1F *	fNhitHisto
	Histogram of number of hits in spill.
TH1F *	fDminHisto0
	Closest approach for particles leaving no hits.
TH1F *	fDminHisto1
	Closest approach for particles leaving hits.
TH2F *	fPhotonAngles
	Photon rate vs angle.
TH2F *	fPhotonAnglesLo
	Photon rate vs angle, low momenta.
TH2F *	fPhotonAnglesMi
	Photon rate vs angle, middle momenta.
TH2F *	fPhotonAnglesHi
	Photon rate vs angle, high momenta.
TH1F *	fPhotonCosQ
	Photon rate vs cos(Q).
TH1F *	fPhotonEnergy
	Photon energy (GeV).
TH2F *	fElectronAngles
	Electron rate vs angle.
TH2F *	fElectronAnglesLo
	Electron rate vs angle, low momenta.
TH2F *	fElectronAnglesMi
	Electron rate vs angle, middle momenta.
TH2F *	fElectronAnglesHi
	Electron rate vs angle, high momenta.
TH1F *	fElectronCosQ
	Electron rate vs cos(Q).
TH1F *	fElectronEnergy
	Electron energy (GeV).
TH2F *	fMuonAngles
	Muon rate vs angle.
TH2F *	fMuonAnglesLo
	Muon rate vs angle, low momenta.
TH2F *	fMuonAnglesMi
	Muon rate vs angle, middle momenta.
TH2F *	fMuonAnglesHi
	Muon rate vs angle, high momenta.
TH1F *	fMuonCosQ
	Muon rate vs cos(Q).
TH1F *	fMuonEnergy
	Muon energy (GeV).
TH2F *	fHitsYZ
	hit locations in YZ view
TH2F *	fHitsXZ
	hit locations in XZ view
CRINFO	fCRInfo
	cosmic ray muon info
TTree *	fTree
	tree holding CRINFO

mcchk::CosmicAna Class Reference

Public Member Functions

Private Attributes

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation