Reference Manual

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