CMSSW Reference Manual

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