Reference Manual

Go to the documentation of this file.
00001 /*
00002  *  See header file for a description of this class.
00003  *
00004  *  $Date: 2008/10/21 10:52:20 $
00005  *  $Revision: 1.2 $
00006  *  \author S. Bolognesi and G. Cerminara - INFN Torino
00007  */
00008 
00009 #include "DTSegment4DQuality.h"
00010 #include "Validation/DTRecHits/interface/DTHitQualityUtils.h"
00011 
00012 #include "DataFormats/MuonDetId/interface/DTWireId.h"
00013 #include "Geometry/DTGeometry/interface/DTChamber.h"
00014 #include "Geometry/DTGeometry/interface/DTGeometry.h"
00015 #include "Geometry/Records/interface/MuonGeometryRecord.h"
00016 
00017 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
00018 #include "DataFormats/DTRecHit/interface/DTRecHitCollection.h"
00019 #include "DataFormats/DTRecHit/interface/DTRecSegment4DCollection.h"
00020 
00021 #include "FWCore/Framework/interface/MakerMacros.h"
00022 #include "FWCore/Framework/interface/Frameworkfwd.h"
00023 #include "FWCore/Framework/interface/Event.h"
00024 #include "FWCore/Framework/interface/ESHandle.h"
00025 
00026 #include "Histograms.h"
00027 
00028 #include "TFile.h"
00029 
00030 #include <iostream>
00031 #include <map>
00032 
00033 
00034 using namespace std;
00035 using namespace edm;
00036 
00037 // Book the histos
00038 
00039 // Constructor
00040 DTSegment4DQuality::DTSegment4DQuality(const ParameterSet& pset)  {
00041   // Get the debug parameter for verbose output
00042   debug = pset.getUntrackedParameter<bool>("debug");
00043   DTHitQualityUtils::debug = debug;
00044 
00045   rootFileName = pset.getUntrackedParameter<string>("rootFileName");
00046   // the name of the simhit collection
00047   simHitLabel = pset.getUntrackedParameter<string>("simHitLabel", "SimG4Object");
00048   // the name of the 4D rec hit collection
00049   segment4DLabel = pset.getUntrackedParameter<string>("segment4DLabel");
00050 
00051   //sigma resolution on position
00052   sigmaResX = pset.getParameter<double>("sigmaResX");
00053   sigmaResY = pset.getParameter<double>("sigmaResY");
00054   //sigma resolution on angle
00055   sigmaResAlpha = pset.getParameter<double>("sigmaResAlpha");
00056   sigmaResBeta = pset.getParameter<double>("sigmaResBeta");
00057 
00058   // Create the root file
00059   theFile = new TFile(rootFileName.c_str(), "RECREATE");
00060   theFile->cd();
00061 
00062   h4DHit= new HRes4DHit ("All");
00063   h4DHit_W0= new HRes4DHit ("W0");
00064   h4DHit_W1= new HRes4DHit ("W1");
00065   h4DHit_W2= new HRes4DHit ("W2");
00066 
00067   hEff_All= new HEff4DHit ("All");
00068   hEff_W0= new HEff4DHit ("W0");
00069   hEff_W1= new HEff4DHit ("W1");
00070   hEff_W2= new HEff4DHit ("W2");
00071 }
00072 
00073 // Destructor
00074 DTSegment4DQuality::~DTSegment4DQuality(){
00075 
00076 }
00077 
00078 void DTSegment4DQuality::endJob() {
00079   // Write the histos to file
00080   theFile->cd();
00081 
00082   h4DHit->Write();
00083   h4DHit_W0->Write();
00084   h4DHit_W1->Write();
00085   h4DHit_W2->Write();
00086 
00087   hEff_All->ComputeEfficiency();
00088   hEff_W0->ComputeEfficiency();
00089   hEff_W1->ComputeEfficiency();
00090   hEff_W2->ComputeEfficiency();
00091 
00092   hEff_All->Write();
00093   hEff_W0->Write();
00094   hEff_W1->Write();
00095   hEff_W2->Write();
00096 
00097   theFile->Close();
00098 } 
00099 
00100 // The real analysis
00101   void DTSegment4DQuality::analyze(const Event & event, const EventSetup& eventSetup){
00102     theFile->cd();
00103 
00104     // Get the DT Geometry
00105     ESHandle<DTGeometry> dtGeom;
00106     eventSetup.get<MuonGeometryRecord>().get(dtGeom);
00107 
00108     // Get the SimHit collection from the event
00109     edm::Handle<PSimHitContainer> simHits;
00110     event.getByLabel(simHitLabel, "MuonDTHits", simHits); //FIXME: second string to be removed
00111 
00112     //Map simHits by chamber
00113     map<DTChamberId, PSimHitContainer > simHitsPerCh;
00114     for(PSimHitContainer::const_iterator simHit = simHits->begin();
00115         simHit != simHits->end(); simHit++){
00116       // Create the id of the chamber (the simHits in the DT known their wireId)
00117       DTChamberId chamberId = (((DTWireId(simHit->detUnitId())).layerId()).superlayerId()).chamberId();
00118       // Fill the map
00119       simHitsPerCh[chamberId].push_back(*simHit);
00120     }
00121 
00122     // Get the 4D rechits from the event
00123     Handle<DTRecSegment4DCollection> segment4Ds;
00124     event.getByLabel(segment4DLabel, segment4Ds);
00125 
00126     // Loop over all chambers containing a segment
00127     DTRecSegment4DCollection::id_iterator chamberId;
00128     for (chamberId = segment4Ds->id_begin();
00129          chamberId != segment4Ds->id_end();
00130          ++chamberId){
00131 
00132       if((*chamberId).station() == 4)
00133         continue; //use DTSegment2DQuality to analyze MB4 performaces
00134 
00135       //------------------------- simHits ---------------------------//
00136       //Get simHits of each chamber
00137       PSimHitContainer simHits =  simHitsPerCh[(*chamberId)];
00138 
00139       // Map simhits per wire
00140       map<DTWireId, PSimHitContainer > simHitsPerWire = DTHitQualityUtils::mapSimHitsPerWire(simHits);
00141       map<DTWireId, const PSimHit*> muSimHitPerWire = DTHitQualityUtils::mapMuSimHitsPerWire(simHitsPerWire);
00142       int nMuSimHit = muSimHitPerWire.size();
00143       if(nMuSimHit == 0 || nMuSimHit == 1) {
00144         if(debug && nMuSimHit == 1)
00145           cout << "[DTSegment4DQuality] Only " << nMuSimHit << " mu SimHit in this chamber, skipping!" << endl;
00146         continue; // If no or only one mu SimHit is found skip this chamber
00147       } 
00148       if(debug)
00149         cout << "=== Chamber " << (*chamberId) << " has " << nMuSimHit << " SimHits" << endl;
00150 
00151       //Find outer and inner mu SimHit to build a segment
00152       pair<const PSimHit*, const PSimHit*> inAndOutSimHit = DTHitQualityUtils::findMuSimSegment(muSimHitPerWire); 
00153 
00154       //Find direction and position of the sim Segment in Chamber RF
00155       pair<LocalVector, LocalPoint> dirAndPosSimSegm = DTHitQualityUtils::findMuSimSegmentDirAndPos(inAndOutSimHit,
00156                                                                                                     (*chamberId),&(*dtGeom));
00157 
00158       LocalVector simSegmLocalDir = dirAndPosSimSegm.first;
00159       LocalPoint simSegmLocalPos = dirAndPosSimSegm.second;
00160       const DTChamber* chamber = dtGeom->chamber(*chamberId);
00161       GlobalPoint simSegmGlobalPos = chamber->toGlobal(simSegmLocalPos);
00162       GlobalVector simSegmGlobalDir = chamber->toGlobal(simSegmLocalDir);
00163 
00164       //phi and theta angle of simulated segment in Chamber RF
00165       float alphaSimSeg = DTHitQualityUtils::findSegmentAlphaAndBeta(simSegmLocalDir).first;
00166       float betaSimSeg = DTHitQualityUtils::findSegmentAlphaAndBeta(simSegmLocalDir).second;
00167       //x,y position of simulated segment in Chamber RF
00168       float xSimSeg = simSegmLocalPos.x(); 
00169       float ySimSeg = simSegmLocalPos.y(); 
00170       //Position (in eta,phi coordinates) in lobal RF
00171       float etaSimSeg = simSegmGlobalPos.eta(); 
00172       float phiSimSeg = simSegmGlobalPos.phi();
00173 
00174       if(debug)
00175         cout<<"  Simulated segment:  local direction "<<simSegmLocalDir<<endl
00176           <<"                      local position  "<<simSegmLocalPos<<endl
00177           <<"                      alpha           "<<alphaSimSeg<<endl
00178           <<"                      beta            "<<betaSimSeg<<endl;
00179 
00180       //---------------------------- recHits --------------------------//
00181       // Get the range of rechit for the corresponding chamberId
00182       bool recHitFound = false;
00183       DTRecSegment4DCollection::range range = segment4Ds->get(*chamberId);
00184       int nsegm = distance(range.first, range.second);
00185       if(debug)
00186         cout << "   Chamber: " << *chamberId << " has " << nsegm
00187           << " 4D segments" << endl;
00188 
00189       if (nsegm!=0) {
00190         // Find the best RecHit: look for the 4D RecHit with the phi angle closest
00191         // to that of segment made of SimHits. 
00192         // RecHits must have delta alpha and delta position within 5 sigma of
00193         // the residual distribution (we are looking for residuals of segments
00194         // usefull to the track fit) for efficency purpose
00195         const DTRecSegment4D* bestRecHit = 0;
00196         bool bestRecHitFound = false;
00197         double deltaAlpha = 99999;
00198 
00199         // Loop over the recHits of this chamberId
00200         for (DTRecSegment4DCollection::const_iterator segment4D = range.first;
00201              segment4D!=range.second;
00202              ++segment4D){
00203           // Check the dimension
00204           if((*segment4D).dimension() != 4) {
00205             if(debug)cout << "[DTSegment4DQuality]***Error: This is not 4D segment!!!" << endl;
00206             continue;
00207           }
00208           // Segment Local Direction and position (in Chamber RF)
00209           LocalVector recSegDirection = (*segment4D).localDirection();
00210           //LocalPoint recSegPosition = (*segment4D).localPosition();
00211 
00212           float recSegAlpha = DTHitQualityUtils::findSegmentAlphaAndBeta(recSegDirection).first;
00213           if(debug)
00214             cout << "  RecSegment direction: " << recSegDirection << endl
00215               << "             position : " <<  (*segment4D).localPosition() << endl
00216               << "             alpha    : " << recSegAlpha << endl
00217               << "             beta     : " <<  DTHitQualityUtils::findSegmentAlphaAndBeta(recSegDirection).second << endl;
00218 
00219           if(fabs(recSegAlpha - alphaSimSeg) < deltaAlpha) {
00220             deltaAlpha = fabs(recSegAlpha - alphaSimSeg);
00221             bestRecHit = &(*segment4D);
00222             bestRecHitFound = true;
00223           }
00224         }  // End of Loop over all 4D RecHits
00225 
00226         if(bestRecHitFound) {
00227           // Best rechit direction and position in Chamber RF
00228           LocalPoint bestRecHitLocalPos = bestRecHit->localPosition();
00229           LocalVector bestRecHitLocalDir = bestRecHit->localDirection();
00230           // Errors on x and y
00231           LocalError bestRecHitLocalPosErr = bestRecHit->localPositionError();
00232           LocalError bestRecHitLocalDirErr = bestRecHit->localDirectionError();
00233 
00234           float alphaBestRHit = DTHitQualityUtils::findSegmentAlphaAndBeta(bestRecHitLocalDir).first;
00235           float betaBestRHit = DTHitQualityUtils::findSegmentAlphaAndBeta(bestRecHitLocalDir).second;
00236           // Errors on alpha and beta
00237 
00238           // Get position and direction using the rx projection (so in SL
00239           // reference frame). Note that x (and y) are swapped wrt to Chamber
00240           // frame
00241           //if (bestRecHit->hasZed() ) {
00242           const DTSLRecSegment2D * zedRecSeg = bestRecHit->zSegment();
00243           LocalPoint bestRecHitLocalPosRZ = zedRecSeg->localPosition();
00244           LocalVector bestRecHitLocalDirRZ = zedRecSeg->localDirection();
00245           // Errors on x and y
00246           LocalError bestRecHitLocalPosErrRZ = zedRecSeg->localPositionError();
00247           LocalError bestRecHitLocalDirErrRZ = zedRecSeg->localDirectionError();
00248 
00249           float alphaBestRHitRZ = DTHitQualityUtils::findSegmentAlphaAndBeta(bestRecHitLocalDirRZ).first;
00250 
00251           // Get SimSeg position and Direction in rZ SL frame 
00252           const DTSuperLayer* sl = dtGeom->superLayer(zedRecSeg->superLayerId());
00253           LocalPoint simSegLocalPosRZTmp = sl->toLocal(simSegmGlobalPos);
00254           LocalVector simSegLocalDirRZ = sl->toLocal(simSegmGlobalDir);
00255           LocalPoint simSegLocalPosRZ =
00256             simSegLocalPosRZTmp + simSegLocalDirRZ*(-simSegLocalPosRZTmp.z()/(cos(simSegLocalDirRZ.theta())));
00257           float alphaSimSegRZ = DTHitQualityUtils::findSegmentAlphaAndBeta(simSegLocalDirRZ).first;
00258 
00259           if (debug) cout <<
00260             "RZ SL: recPos " << bestRecHitLocalPosRZ <<
00261               "recDir " << bestRecHitLocalDirRZ <<
00262               "recAlpha " << alphaBestRHitRZ << endl <<
00263               "RZ SL: simPos " << simSegLocalPosRZ <<
00264               "simDir " << simSegLocalDirRZ <<
00265               "simAlpha " << alphaSimSegRZ << endl ;
00266           //}
00267 
00268 
00269           if(fabs(alphaBestRHit - alphaSimSeg) < 5*sigmaResAlpha &&
00270              fabs(betaBestRHit - betaSimSeg) < 5*sigmaResBeta &&
00271              fabs(bestRecHitLocalPos.x() - xSimSeg) < 5*sigmaResX &&
00272              fabs(bestRecHitLocalPos.y() - ySimSeg) < 5*sigmaResY) {
00273             recHitFound = true;
00274           }
00275 
00276           // Fill Residual histos
00277           HRes4DHit *histo=0;
00278 
00279           if((*chamberId).wheel() == 0)
00280             histo = h4DHit_W0;
00281           else if(abs((*chamberId).wheel()) == 1)
00282             histo = h4DHit_W1;
00283           else if(abs((*chamberId).wheel()) == 2)
00284             histo = h4DHit_W2;
00285 
00286           histo->Fill(alphaSimSeg,
00287                       alphaBestRHit,
00288                       betaSimSeg,
00289                       betaBestRHit,
00290                       xSimSeg, 
00291                       bestRecHitLocalPos.x(),
00292                       ySimSeg,
00293                       bestRecHitLocalPos.y(),
00294                       etaSimSeg, 
00295                       phiSimSeg,
00296                       bestRecHitLocalPosRZ.x(),
00297                       simSegLocalPosRZ.x(),
00298                       alphaBestRHitRZ,
00299                       alphaSimSegRZ,
00300                       sqrt(bestRecHitLocalDirErr.xx()),
00301                       sqrt(bestRecHitLocalDirErr.yy()),
00302                       sqrt(bestRecHitLocalPosErr.xx()),
00303                       sqrt(bestRecHitLocalPosErr.yy()),
00304                       sqrt(bestRecHitLocalDirErrRZ.xx()),
00305                       sqrt(bestRecHitLocalPosErrRZ.xx())
00306                      );
00307 
00308           h4DHit->Fill(alphaSimSeg,
00309                        alphaBestRHit,
00310                        betaSimSeg,
00311                        betaBestRHit,
00312                        xSimSeg, 
00313                        bestRecHitLocalPos.x(),
00314                        ySimSeg,
00315                        bestRecHitLocalPos.y(),
00316                        etaSimSeg, 
00317                        phiSimSeg,
00318                        bestRecHitLocalPosRZ.x(),
00319                        simSegLocalPosRZ.x(),
00320                        alphaBestRHitRZ,
00321                        alphaSimSegRZ,
00322                        sqrt(bestRecHitLocalDirErr.xx()),
00323                        sqrt(bestRecHitLocalDirErr.yy()),
00324                        sqrt(bestRecHitLocalPosErr.xx()),
00325                        sqrt(bestRecHitLocalPosErr.yy()),
00326                        sqrt(bestRecHitLocalDirErrRZ.xx()),
00327                        sqrt(bestRecHitLocalPosErrRZ.xx())
00328                       );
00329         }
00330       } //end of if(nsegm!=0)
00331 
00332       // Fill Efficiency plot
00333       HEff4DHit *heff = 0;
00334 
00335       if((*chamberId).wheel() == 0)
00336         heff = hEff_W0;
00337       else if(abs((*chamberId).wheel()) == 1)
00338         heff = hEff_W1;
00339       else if(abs((*chamberId).wheel()) == 2)
00340         heff = hEff_W2;
00341       heff->Fill(etaSimSeg, phiSimSeg, xSimSeg, ySimSeg, alphaSimSeg, betaSimSeg, recHitFound);
00342       hEff_All->Fill(etaSimSeg, phiSimSeg, xSimSeg, ySimSeg, alphaSimSeg, betaSimSeg, recHitFound);
00343     } // End of loop over chambers
00344   }
/data/refman/pasoursint/CMSSW_4_1_8_patch13/src/Validation/DTRecHits/plugins/DTSegment4DQuality.cc
