Reference Manual

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