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.8 $
00006  *  \author S. Bolognesi and G. Cerminara - INFN Torino
00007  */
00008 
00009 #include "DTSegment2DQuality.h"
00010 #include "Validation/DTRecHits/interface/DTHitQualityUtils.h"
00011 
00012 #include "Geometry/DTGeometry/interface/DTSuperLayer.h"
00013 #include "DataFormats/MuonDetId/interface/DTWireId.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/DTRecSegment2DCollection.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 "DQMServices/Core/interface/DQMStore.h"
00031 #include "DQMServices/Core/interface/MonitorElement.h"
00032 #include <iostream>
00033 #include <map>
00034 
00035 using namespace std;
00036 using namespace edm;
00037 
00038 // Constructor
00039 DTSegment2DQuality::DTSegment2DQuality(const ParameterSet& pset)  {
00040    // Get the debug parameter for verbose output
00041   debug = pset.getUntrackedParameter<bool>("debug");
00042   DTHitQualityUtils::debug = debug;
00043   // the name of the simhit collection
00044   simHitLabel = pset.getUntrackedParameter<InputTag>("simHitLabel");
00045   // the name of the 2D rec hit collection
00046   segment2DLabel = pset.getUntrackedParameter<InputTag>("segment2DLabel");
00047 
00048   //sigma resolution on position
00049   sigmaResPos = pset.getParameter<double>("sigmaResPos");
00050   //sigma resolution on angle
00051   sigmaResAngle = pset.getParameter<double>("sigmaResAngle");
00052 
00053   // Create the root file
00054   //theFile = new TFile(rootFileName.c_str(), "RECREATE");
00055   //theFile->cd();
00056 
00057   if(debug)
00058     cout << "[DTSegment2DQuality] Constructor called " << endl;
00059   
00060   // ----------------------                 
00061   // get hold of back-end interface 
00062   dbe_ = 0;
00063   dbe_ = Service<DQMStore>().operator->();
00064   if ( dbe_ ) {
00065     if (debug) {
00066       dbe_->setVerbose(1);
00067     } else {
00068       dbe_->setVerbose(0);
00069     }
00070   }
00071   if ( dbe_ ) {
00072     if ( debug ) dbe_->showDirStructure();
00073   }
00074 
00075   h2DHitRPhi = new HRes2DHit ("RPhi",dbe_,true,true);
00076   h2DHitRZ= new HRes2DHit ("RZ",dbe_,true,true);
00077   h2DHitRZ_W0= new HRes2DHit ("RZ_W0",dbe_,true,true);
00078   h2DHitRZ_W1= new HRes2DHit ("RZ_W1",dbe_,true,true);
00079   h2DHitRZ_W2= new HRes2DHit ("RZ_W2",dbe_,true,true);
00080 
00081   h2DHitEff_RPhi= new HEff2DHit ("RPhi",dbe_);
00082   h2DHitEff_RZ= new HEff2DHit ("RZ",dbe_);
00083   h2DHitEff_RZ_W0= new HEff2DHit ("RZ_W0",dbe_);
00084   h2DHitEff_RZ_W1= new HEff2DHit ("RZ_W1",dbe_);
00085   h2DHitEff_RZ_W2= new HEff2DHit ("RZ_W2",dbe_);
00086   if(debug)
00087     cout << "[DTSegment2DQuality] hitsos created " << endl;
00088 }
00089 
00090 // Destructor
00091 DTSegment2DQuality::~DTSegment2DQuality(){
00092 
00093 }
00094 
00095 void DTSegment2DQuality::endJob() {
00096   // Write the histos to file
00097   //theFile->cd();
00098 
00099   //h2DHitRPhi->Write();
00100   //h2DHitRZ->Write();
00101   //h2DHitRZ_W0->Write();
00102   //h2DHitRZ_W1->Write();
00103   //h2DHitRZ_W2->Write();
00104 
00105   h2DHitEff_RPhi->ComputeEfficiency();
00106   h2DHitEff_RZ->ComputeEfficiency();
00107   h2DHitEff_RZ_W0->ComputeEfficiency();
00108   h2DHitEff_RZ_W1->ComputeEfficiency();
00109   h2DHitEff_RZ_W2->ComputeEfficiency();
00110 
00111   //h2DHitEff_RPhi->Write();
00112   //h2DHitEff_RZ->Write();
00113   //h2DHitEff_RZ_W0->Write();
00114   //h2DHitEff_RZ_W1->Write();
00115   //h2DHitEff_RZ_W2->Write();
00116   //if ( rootFileName.size() != 0 && dbe_ ) dbe_->save(rootFileName); 
00117 
00118   //theFile->Close();
00119 } 
00120 
00121 // The real analysis
00122 void DTSegment2DQuality::analyze(const Event & event, const EventSetup& eventSetup){
00123   //theFile->cd();
00124 
00125   // Get the DT Geometry
00126   ESHandle<DTGeometry> dtGeom;
00127   eventSetup.get<MuonGeometryRecord>().get(dtGeom);
00128 
00129   // Get the SimHit collection from the event
00130   edm::Handle<PSimHitContainer> simHits;
00131   event.getByLabel(simHitLabel, simHits); //FIXME: second string to be removed
00132 
00133   //Map simHits by sl
00134   map<DTSuperLayerId, PSimHitContainer > simHitsPerSl;
00135   for(PSimHitContainer::const_iterator simHit = simHits->begin();
00136       simHit != simHits->end(); simHit++){
00137     // Create the id of the sl (the simHits in the DT known their wireId)
00138     DTSuperLayerId slId = ((DTWireId(simHit->detUnitId())).layerId()).superlayerId();
00139     // Fill the map
00140     simHitsPerSl[slId].push_back(*simHit);
00141   }
00142 
00143   // Get the 2D rechits from the event
00144   Handle<DTRecSegment2DCollection> segment2Ds;
00145   event.getByLabel(segment2DLabel, segment2Ds);
00146 
00147   if(!segment2Ds.isValid()) {
00148     if(debug) cout << "[DTSegment2DQuality]**Warning: no 2DSegments with label: " << segment2DLabel
00149                    << " in this event, skipping!" << endl;
00150     return;
00151   }    
00152 
00153   // Loop over all superlayers containing a segment
00154   DTRecSegment2DCollection::id_iterator slId;
00155   for (slId = segment2Ds->id_begin();
00156        slId != segment2Ds->id_end();
00157        ++slId){
00158 
00159       //------------------------- simHits ---------------------------//
00160     //Get simHits of each superlayer
00161     PSimHitContainer simHits =  simHitsPerSl[(*slId)];
00162        
00163     // Map simhits per wire
00164     map<DTWireId, PSimHitContainer > simHitsPerWire = DTHitQualityUtils::mapSimHitsPerWire(simHits);
00165     map<DTWireId, const PSimHit*> muSimHitPerWire = DTHitQualityUtils::mapMuSimHitsPerWire(simHitsPerWire);
00166     int nMuSimHit = muSimHitPerWire.size();
00167     if(nMuSimHit == 0 || nMuSimHit == 1) {
00168       if(debug && nMuSimHit == 1)
00169         cout << "[DTSegment2DQuality] Only " << nMuSimHit << " mu SimHit in this SL, skipping!" << endl;
00170       continue; // If no or only one mu SimHit is found skip this SL
00171     } 
00172     if(debug)
00173       cout << "=== SL " << (*slId) << " has " << nMuSimHit << " SimHits" << endl;
00174   
00175     //Find outer and inner mu SimHit to build a segment
00176     pair<const PSimHit*, const PSimHit*> inAndOutSimHit = DTHitQualityUtils::findMuSimSegment(muSimHitPerWire); 
00177     //Check that outermost and innermost SimHit are not the same
00178     if(inAndOutSimHit.first ==inAndOutSimHit.second ) {
00179       cout << "[DTHitQualityUtils]***Warning: outermost and innermost SimHit are the same!" << endl;
00180       continue;
00181     }
00182 
00183     //Find direction and position of the sim Segment in SL RF
00184     pair<LocalVector, LocalPoint> dirAndPosSimSegm = DTHitQualityUtils::findMuSimSegmentDirAndPos(inAndOutSimHit,
00185                                                                                                   (*slId),&(*dtGeom));
00186 
00187     LocalVector simSegmLocalDir = dirAndPosSimSegm.first;
00188     LocalPoint simSegmLocalPos = dirAndPosSimSegm.second;
00189     if(debug)
00190       cout<<"  Simulated segment:  local direction "<<simSegmLocalDir<<endl
00191         <<"                      local position  "<<simSegmLocalPos<<endl;
00192     const DTSuperLayer* superLayer = dtGeom->superLayer(*slId);
00193     GlobalPoint simSegmGlobalPos = superLayer->toGlobal(simSegmLocalPos);
00194 
00195     //Atan(x/z) angle and x position in SL RF
00196     float angleSimSeg = DTHitQualityUtils::findSegmentAlphaAndBeta(simSegmLocalDir).first;
00197     float posSimSeg = simSegmLocalPos.x(); 
00198     //Position (in eta,phi coordinates) in the global RF
00199     float etaSimSeg = simSegmGlobalPos.eta(); 
00200     float phiSimSeg = simSegmGlobalPos.phi();
00201 
00202 
00203     //---------------------------- recHits --------------------------//
00204     // Get the range of rechit for the corresponding slId
00205     bool recHitFound = false;
00206     DTRecSegment2DCollection::range range = segment2Ds->get(*slId);
00207     int nsegm = distance(range.first, range.second);
00208     if(debug)
00209       cout << "   Sl: " << *slId << " has " << nsegm
00210            << " 2D segments" << endl;
00211 
00212     if (nsegm!=0) {
00213       // Find the best RecHit: look for the 2D RecHit with the angle closest
00214       // to that of segment made of SimHits. 
00215       // RecHits must have delta alpha and delta position within 5 sigma of
00216       // the residual distribution (we are looking for residuals of segments
00217       // usefull to the track fit) for efficency purpose
00218       const DTRecSegment2D* bestRecHit = 0;
00219       bool bestRecHitFound = false;
00220       double deltaAlpha = 99999;
00221   
00222       // Loop over the recHits of this slId
00223       for (DTRecSegment2DCollection::const_iterator segment2D = range.first;
00224            segment2D!=range.second;
00225            ++segment2D){
00226         // Check the dimension
00227         if((*segment2D).dimension() != 2) {
00228       if(debug) cout << "[DTSegment2DQuality]***Error: This is not 2D segment!!!" << endl;
00229           abort();
00230         }
00231         // Segment Local Direction and position (in SL RF)
00232         LocalVector recSegDirection = (*segment2D).localDirection();
00233         LocalPoint recSegPosition = (*segment2D).localPosition();
00234       
00235         float recSegAlpha = DTHitQualityUtils::findSegmentAlphaAndBeta(recSegDirection).first;
00236         if(debug)
00237       cout << "  RecSegment direction: " << recSegDirection << endl
00238              << "             position : " << recSegPosition << endl
00239              << "             alpha    : " << recSegAlpha << endl;
00240       
00241         if(fabs(recSegAlpha - angleSimSeg) < deltaAlpha) {
00242           deltaAlpha = fabs(recSegAlpha - angleSimSeg);
00243           bestRecHit = &(*segment2D);
00244           bestRecHitFound = true;
00245         }
00246       }  // End of Loop over all 2D RecHits
00247 
00248       if(bestRecHitFound) {
00249         // Best rechit direction and position in SL RF
00250         LocalPoint bestRecHitLocalPos = bestRecHit->localPosition();
00251         LocalVector bestRecHitLocalDir = bestRecHit->localDirection();
00252 
00253     LocalError bestRecHitLocalPosErr = bestRecHit->localPositionError();
00254     LocalError bestRecHitLocalDirErr = bestRecHit->localDirectionError();
00255 
00256         float angleBestRHit = DTHitQualityUtils::findSegmentAlphaAndBeta(bestRecHitLocalDir).first;
00257 
00258         if(fabs(angleBestRHit - angleSimSeg) < 5*sigmaResAngle &&
00259            fabs(bestRecHitLocalPos.x() - posSimSeg) < 5*sigmaResPos) {
00260           recHitFound = true;
00261         }
00262 
00263         // Fill Residual histos
00264         HRes2DHit *hRes = 0;
00265         if((*slId).superlayer() == 1 || (*slId).superlayer() == 3) { //RPhi SL
00266           hRes = h2DHitRPhi;
00267         } else if((*slId).superlayer() == 2) { //RZ SL
00268           h2DHitRZ->Fill(angleSimSeg,
00269                     angleBestRHit,
00270                     posSimSeg,
00271                     bestRecHitLocalPos.x(),
00272                     etaSimSeg,
00273                     phiSimSeg,
00274                     sqrt(bestRecHitLocalPosErr.xx()),
00275                     sqrt(bestRecHitLocalDirErr.xx())
00276                     );
00277           if(abs((*slId).wheel()) == 0)
00278             hRes = h2DHitRZ_W0;
00279           else if(abs((*slId).wheel()) == 1)
00280             hRes = h2DHitRZ_W1;
00281           else if(abs((*slId).wheel()) == 2)
00282             hRes = h2DHitRZ_W2;
00283         }
00284         hRes->Fill(angleSimSeg,
00285                angleBestRHit,
00286                posSimSeg,
00287                bestRecHitLocalPos.x(),
00288                etaSimSeg,
00289                phiSimSeg,
00290                sqrt(bestRecHitLocalPosErr.xx()),
00291                sqrt(bestRecHitLocalDirErr.xx())
00292                );
00293       }
00294     } //end of if(nsegm!=0)
00295 
00296       // Fill Efficiency plot
00297     HEff2DHit *hEff = 0;
00298     if((*slId).superlayer() == 1 || (*slId).superlayer() == 3) { //RPhi SL
00299       hEff = h2DHitEff_RPhi;
00300     } else if((*slId).superlayer() == 2) { //RZ SL
00301       h2DHitEff_RZ->Fill(etaSimSeg, phiSimSeg, posSimSeg, angleSimSeg, recHitFound);
00302       if(abs((*slId).wheel()) == 0)
00303         hEff = h2DHitEff_RZ_W0;
00304       else if(abs((*slId).wheel()) == 1)
00305         hEff = h2DHitEff_RZ_W1;
00306       else if(abs((*slId).wheel()) == 2)
00307         hEff = h2DHitEff_RZ_W2;
00308     }
00309     hEff->Fill(etaSimSeg, phiSimSeg, posSimSeg, angleSimSeg, recHitFound);
00310   } // End of loop over superlayers
00311 }
00312 
00313 
00314 
00315 
00316 
00317
/data/refman/pasoursint/CMSSW_4_4_5_patch3/src/Validation/DTRecHits/plugins/DTSegment2DQuality.cc
