Reference Manual

Go to the documentation of this file.
00001 
00002 /*
00003  *  DQM muon alignment analysis monitoring
00004  *
00005  *  $Date: 2010/03/29 13:18:44 $
00006  *  $Revision: 1.6 $
00007  *  \author J. Fernandez - Univ. Oviedo <Javier.Fernandez@cern.ch>
00008  */
00009 
00010 #include "DQMOffline/Alignment/interface/MuonAlignment.h"
00011 
00012 
00013 MuonAlignment::MuonAlignment(const edm::ParameterSet& pSet) {
00014 
00015     metname = "MuonAlignment";
00016 
00017     LogTrace(metname)<<"[MuonAlignment] Constructor called!"<<std::endl;
00018 
00019     parameters = pSet;
00020 
00021     theMuonCollectionLabel = parameters.getParameter<edm::InputTag>("MuonCollection");
00022   
00023     theRecHits4DTagDT = parameters.getParameter<edm::InputTag>("RecHits4DDTCollectionTag");
00024     theRecHits4DTagCSC = parameters.getParameter<edm::InputTag>("RecHits4DCSCCollectionTag");
00025   
00026     resLocalXRangeStation1 = parameters.getUntrackedParameter<double>("resLocalXRangeStation1");
00027     resLocalXRangeStation2 = parameters.getUntrackedParameter<double>("resLocalXRangeStation2");
00028     resLocalXRangeStation3 = parameters.getUntrackedParameter<double>("resLocalXRangeStation3");
00029     resLocalXRangeStation4 = parameters.getUntrackedParameter<double>("resLocalXRangeStation4");
00030     resLocalYRangeStation1 = parameters.getUntrackedParameter<double>("resLocalYRangeStation1");
00031     resLocalYRangeStation2 = parameters.getUntrackedParameter<double>("resLocalYRangeStation2");
00032     resLocalYRangeStation3 = parameters.getUntrackedParameter<double>("resLocalYRangeStation3");
00033     resLocalYRangeStation4 = parameters.getUntrackedParameter<double>("resLocalYRangeStation4");
00034     resPhiRange = parameters.getUntrackedParameter<double>("resPhiRange");
00035     resThetaRange = parameters.getUntrackedParameter<double>("resThetaRange");
00036 
00037     meanPositionRange = parameters.getUntrackedParameter<double>("meanPositionRange");
00038     rmsPositionRange = parameters.getUntrackedParameter<double>("rmsPositionRange");
00039     meanAngleRange = parameters.getUntrackedParameter<double>("meanAngleRange");
00040     rmsAngleRange = parameters.getUntrackedParameter<double>("rmsAngleRange");
00041 
00042     nbins = parameters.getUntrackedParameter<unsigned int>("nbins");
00043     min1DTrackRecHitSize = parameters.getUntrackedParameter<unsigned int>("min1DTrackRecHitSize");
00044     min4DTrackSegmentSize = parameters.getUntrackedParameter<unsigned int>("min4DTrackSegmentSize");
00045 
00046     doDT = parameters.getUntrackedParameter<bool>("doDT");
00047     doCSC = parameters.getUntrackedParameter<bool>("doCSC");
00048     doSummary = parameters.getUntrackedParameter<bool>("doSummary");
00049 
00050     numberOfTracks=0;
00051     numberOfHits=0;
00052   
00053     MEFolderName = parameters.getParameter<std::string>("FolderName");  
00054     topFolder << MEFolderName+"/Alignment/Muon";
00055 }
00056 
00057 MuonAlignment::~MuonAlignment() { 
00058 }
00059 
00060 
00061 void MuonAlignment::beginJob() {
00062 
00063 
00064     LogTrace(metname)<<"[MuonAlignment] Parameters initialization";
00065   
00066     if(!(doDT || doCSC) ) { 
00067         edm::LogError("MuonAlignment") <<" Error!! At least one Muon subsystem (DT or CSC) must be monitorized!!" << std::endl;
00068         edm::LogError("MuonAlignment") <<" Please enable doDT or doCSC to True in your python cfg file!!!" << std::endl;
00069         exit(1);
00070     }
00071   
00072     dbe = edm::Service<DQMStore>().operator->();
00073 
00074     if (doSummary){
00075         if (doDT){
00076             dbe->setCurrentFolder(topFolder.str()+"/DT");
00077             hLocalPositionDT=dbe->book2D("hLocalPositionDT","Local DT position (cm) absolute MEAN residuals;Sector;;cm", 14,1, 15,40,0,40);
00078             hLocalAngleDT=dbe->book2D("hLocalAngleDT","Local DT angle (rad) absolute MEAN residuals;Sector;;rad", 14,1, 15,40,0,40); 
00079             hLocalPositionRmsDT=dbe->book2D("hLocalPositionRmsDT","Local DT position (cm) RMS residuals;Sector;;cm", 14,1, 15,40,0,40);
00080             hLocalAngleRmsDT=dbe->book2D("hLocalAngleRmsDT","Local DT angle (rad) RMS residuals;Sector;;rad", 14,1, 15,40,0,40); 
00081 
00082             hLocalXMeanDT=dbe->book1D("hLocalXMeanDT","Distribution of absolute MEAN Local X (cm) residuals for DT;<X> (cm);number of chambers",100,0,meanPositionRange);
00083             hLocalXRmsDT=dbe->book1D("hLocalXRmsDT","Distribution of RMS Local X (cm) residuals for DT;X RMS (cm);number of chambers", 100,0,rmsPositionRange);
00084             hLocalYMeanDT=dbe->book1D("hLocalYMeanDT","Distribution of absolute MEAN Local Y (cm) residuals for DT;<Y> (cm);number of chambers", 100,0,meanPositionRange);
00085             hLocalYRmsDT=dbe->book1D("hLocalYRmsDT","Distribution of RMS Local Y (cm) residuals for DT;Y RMS (cm);number of chambers", 100,0,rmsPositionRange);
00086 
00087             hLocalPhiMeanDT=dbe->book1D("hLocalPhiMeanDT","Distribution of MEAN #phi (rad) residuals for DT;<#phi>(rad);number of chambers", 100,-meanAngleRange,meanAngleRange);
00088             hLocalPhiRmsDT=dbe->book1D("hLocalPhiRmsDT","Distribution of RMS #phi (rad) residuals for DT;#phi RMS (rad);number of chambers", 100,0,rmsAngleRange);
00089             hLocalThetaMeanDT=dbe->book1D("hLocalThetaMeanDT","Distribution of MEAN #theta (rad) residuals for DT;<#theta>(rad);number of chambers", 100,-meanAngleRange,meanAngleRange);
00090             hLocalThetaRmsDT=dbe->book1D("hLocalThetaRmsDT","Distribution of RMS #theta (rad) residuals for DT;#theta RMS (rad);number of chambers",100,0,rmsAngleRange);
00091         }
00092         
00093         if (doCSC){
00094             dbe->setCurrentFolder(topFolder.str()+"/CSC");
00095             hLocalPositionCSC=dbe->book2D("hLocalPositionCSC","Local CSC position (cm) absolute MEAN residuals;Sector;;cm",36,1,37,40,0,40);
00096             hLocalAngleCSC=dbe->book2D("hLocalAngleCSC","Local CSC angle (rad) absolute MEAN residuals;Sector;;rad", 36,1,37,40,0,40); 
00097             hLocalPositionRmsCSC=dbe->book2D("hLocalPositionRmsCSC","Local CSC position (cm) RMS residuals;Sector;;cm", 36,1,37,40,0,40);
00098             hLocalAngleRmsCSC=dbe->book2D("hLocalAngleRmsCSC","Local CSC angle (rad) RMS residuals;Sector;;rad", 36,1,37,40,0,40); 
00099         
00100             hLocalXMeanCSC=dbe->book1D("hLocalXMeanCSC","Distribution of absolute MEAN Local X (cm) residuals for CSC;<X> (cm);number of chambers",100,0,meanPositionRange);
00101             hLocalXRmsCSC=dbe->book1D("hLocalXRmsCSC","Distribution of RMS Local X (cm) residuals for CSC;X RMS (cm);number of chambers", 100,0,rmsPositionRange);
00102             hLocalYMeanCSC=dbe->book1D("hLocalYMeanCSC","Distribution of absolute MEAN Local Y (cm) residuals for CSC;<Y> (cm);number of chambers", 100,0,meanPositionRange);
00103             hLocalYRmsCSC=dbe->book1D("hLocalYRmsCSC","Distribution of RMS Local Y (cm) residuals for CSC;Y RMS (cm);number of chambers", 100,0,rmsPositionRange);
00104 
00105             hLocalPhiMeanCSC=dbe->book1D("hLocalPhiMeanCSC","Distribution of absolute MEAN #phi (rad) residuals for CSC;<#phi>(rad);number of chambers", 100,0,meanAngleRange);
00106             hLocalPhiRmsCSC=dbe->book1D("hLocalPhiRmsCSC","Distribution of RMS #phi (rad) residuals for CSC;#phi RMS (rad);number of chambers", 100,0,rmsAngleRange);
00107             hLocalThetaMeanCSC=dbe->book1D("hLocalThetaMeanCSC","Distribution of absolute MEAN #theta (rad) residuals for CSC;<#theta>(rad);number of chambers", 100,0,meanAngleRange);
00108             hLocalThetaRmsCSC=dbe->book1D("hLocalThetaRmsCSC","Distribution of RMS #theta (rad) residuals for CSC;#theta RMS (rad);number of chambers",100,0,rmsAngleRange);
00109         }
00110     }
00111 
00112 
00113         // Chamber individual histograms
00114         // I need to create all of them even if they are empty to allow proper root merging
00115 
00116         // variables for histos ranges  
00117         double rangeX=0,rangeY=0;
00118         std::string nameOfHistoLocalX,nameOfHistoLocalY,nameOfHistoLocalPhi,nameOfHistoLocalTheta;
00119 
00120         for (int station = -4; station<5; station++){
00121 
00122         //This piece of code calculates the range of the residuals
00123         switch(abs(station)) {
00124             case 1:
00125             {rangeX = resLocalXRangeStation1; rangeY = resLocalYRangeStation1;}
00126             break;
00127             case 2:
00128             {rangeX = resLocalXRangeStation2; rangeY = resLocalYRangeStation2;}
00129             break;
00130             case 3:
00131             {rangeX = resLocalXRangeStation3; rangeY = resLocalYRangeStation3;}
00132             break;
00133             case 4:
00134             {rangeX = resLocalXRangeStation4; rangeY = resLocalYRangeStation4;}
00135             break;
00136             default:
00137                 break;
00138         }
00139         if (doDT){
00140             if(station>0){
00141                                         
00142                 for(int wheel=-2;wheel<3;wheel++){
00143                         
00144                     for (int sector=1;sector<15;sector++){
00145                                 
00146                         if(!((sector==13 || sector ==14) && station!=4)){
00147                                         
00148                             std::stringstream Wheel; Wheel<<wheel;
00149                             std::stringstream Station; Station<<station;
00150                             std::stringstream Sector; Sector<<sector;
00151                                                                                 
00152                             nameOfHistoLocalX="ResidualLocalX_W"+Wheel.str()+"MB"+Station.str()+"S"+Sector.str();
00153                             nameOfHistoLocalPhi= "ResidualLocalPhi_W"+Wheel.str()+"MB"+Station.str()+"S"+Sector.str();
00154                             nameOfHistoLocalTheta= "ResidualLocalTheta_W"+Wheel.str()+"MB"+Station.str()+"S"+Sector.str();
00155                             nameOfHistoLocalY= "ResidualLocalY_W"+Wheel.str()+"MB"+Station.str()+"S"+Sector.str();
00156                                                                                
00157                             dbe->setCurrentFolder(topFolder.str()+
00158                                                   "/DT/Wheel"+Wheel.str()+
00159                                                   "/Station"+Station.str()+
00160                                                   "/Sector"+Sector.str());
00161 
00162                             //Create ME and push histos into their respective vectors
00163 
00164                             MonitorElement *histoLocalX = dbe->book1D(nameOfHistoLocalX, nameOfHistoLocalX, nbins, -rangeX, rangeX);
00165                             unitsLocalX.push_back(histoLocalX);
00166                             MonitorElement *histoLocalPhi = dbe->book1D(nameOfHistoLocalPhi, nameOfHistoLocalPhi, nbins, -resPhiRange, resPhiRange);
00167                             unitsLocalPhi.push_back(histoLocalPhi);
00168                             MonitorElement *histoLocalTheta = dbe->book1D(nameOfHistoLocalTheta, nameOfHistoLocalTheta, nbins, -resThetaRange, resThetaRange);
00169                             unitsLocalTheta.push_back(histoLocalTheta);
00170                             MonitorElement *histoLocalY = dbe->book1D(nameOfHistoLocalY, nameOfHistoLocalY, nbins, -rangeY, rangeY);
00171                             unitsLocalY.push_back(histoLocalY);
00172                                         }
00173                     }
00174                 }
00175             } //station>0
00176         }// doDT
00177         
00178         if (doCSC){
00179             if(station!=0){
00180 
00181                 for(int ring=1;ring<5;ring++){
00182 
00183                     for(int chamber=1;chamber<37;chamber++){
00184                                 
00185                         if( !( ((abs(station)==2 || abs(station)==3 || abs(station)==4) && ring==1 && chamber>18) || 
00186                                ((abs(station)==2 || abs(station)==3 || abs(station)==4) && ring>2)) ){
00187                             std::stringstream Ring; Ring<<ring;
00188                             std::stringstream Station; Station<<station;
00189                             std::stringstream Chamber; Chamber<<chamber;
00190                                                                                
00191                             nameOfHistoLocalX="ResidualLocalX_ME"+Station.str()+"R"+Ring.str()+"C"+Chamber.str();
00192                             nameOfHistoLocalPhi= "ResidualLocalPhi_ME"+Station.str()+"R"+Ring.str()+"C"+Chamber.str();
00193                             nameOfHistoLocalTheta= "ResidualLocalTheta_ME"+Station.str()+"R"+Ring.str()+"C"+Chamber.str();
00194                             nameOfHistoLocalY= "ResidualLocalY_ME"+Station.str()+"R"+Ring.str()+"C"+Chamber.str();
00195                                                                                
00196                             dbe->setCurrentFolder(topFolder.str()+
00197                                                   "/CSC/Station"+Station.str()+
00198                                                   "/Ring"+Ring.str()+
00199                                                   "/Chamber"+Chamber.str());
00200 
00201                             //Create ME and push histos into their respective vectors
00202 
00203                             MonitorElement *histoLocalX = dbe->book1D(nameOfHistoLocalX, nameOfHistoLocalX, nbins, -rangeX, rangeX);
00204                             unitsLocalX.push_back(histoLocalX);
00205                             MonitorElement *histoLocalPhi = dbe->book1D(nameOfHistoLocalPhi, nameOfHistoLocalPhi, nbins, -resPhiRange, resPhiRange);
00206                             unitsLocalPhi.push_back(histoLocalPhi);
00207                             MonitorElement *histoLocalTheta = dbe->book1D(nameOfHistoLocalTheta, nameOfHistoLocalTheta, nbins, -resThetaRange, resThetaRange);
00208                             unitsLocalTheta.push_back(histoLocalTheta);
00209                             MonitorElement *histoLocalY = dbe->book1D(nameOfHistoLocalY, nameOfHistoLocalY, nbins, -rangeY, rangeY);
00210                             unitsLocalY.push_back(histoLocalY);
00211 
00212                         }
00213                     }
00214                 }
00215             } // station!=0
00216         }// doCSC
00217 
00218         } // loop on stations
00219 }
00220 
00221 
00222 void MuonAlignment::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
00223 
00224     LogTrace(metname)<<"[MuonAlignment] Analysis of event # ";
00225   
00226     edm::ESHandle<MagneticField> theMGField;
00227     iSetup.get<IdealMagneticFieldRecord>().get(theMGField);
00228 
00229     edm::ESHandle<GlobalTrackingGeometry> theTrackingGeometry;
00230     iSetup.get<GlobalTrackingGeometryRecord>().get(theTrackingGeometry);
00231 
00232         edm::ESHandle<Propagator> thePropagatorOpp;
00233         iSetup.get<TrackingComponentsRecord>().get( "SmartPropagatorOpposite", thePropagatorOpp );
00234   
00235         edm::ESHandle<Propagator> thePropagatorAlo;
00236         iSetup.get<TrackingComponentsRecord>().get( "SmartPropagator", thePropagatorAlo );
00237 
00238 //      edm::ESHandle<Propagator> thePropagator;
00239 //      iSetup.get<TrackingComponentsRecord>().get( "SmartPropagatorAnyOpposite", thePropagator );
00240 
00241     // Get the RecoMuons collection from the event
00242     edm::Handle<reco::TrackCollection> muons;
00243     iEvent.getByLabel(theMuonCollectionLabel, muons);
00244 
00245     // Get the 4D DTSegments
00246     edm::Handle<DTRecSegment4DCollection> all4DSegmentsDT;
00247     iEvent.getByLabel(theRecHits4DTagDT, all4DSegmentsDT);
00248     DTRecSegment4DCollection::const_iterator segmentDT;
00249 
00250     // Get the 4D CSCSegments
00251     edm::Handle<CSCSegmentCollection> all4DSegmentsCSC;
00252     iEvent.getByLabel(theRecHits4DTagCSC, all4DSegmentsCSC);
00253     CSCSegmentCollection::const_iterator segmentCSC;
00254   
00255     //Vectors used to perform the matching between Segments and hits from Track
00256     intDVector indexCollectionDT;
00257     intDVector indexCollectionCSC;
00258 
00259         
00260 //        thePropagator = new SteppingHelixPropagator(&*theMGField, alongMomentum);
00261 
00262     int countTracks   = 0;
00263     reco::TrackCollection::const_iterator muon;
00264     for (muon = muons->begin(); muon != muons->end(); ++muon){
00265 //      if(muon->isGlobalMuon()){
00266 //      if(muon->isStandAloneMuon()){
00267 
00268         int countPoints   = 0;
00269 
00270 //      reco::TrackRef trackTR = muon->innerTrack();
00271 //      reco::TrackRef trackSA = muon->outerTrack();
00272         //reco::TrackRef trackSA = muon;
00273     
00274         if(muon->recHitsSize()>(min1DTrackRecHitSize-1)) {
00275 
00276 //      reco::TransientTrack tTrackTR( *trackTR, &*theMGField, theTrackingGeometry );
00277             reco::TransientTrack tTrackSA(*muon,&*theMGField,theTrackingGeometry); 
00278 
00279             // Adapted code for muonCosmics
00280    
00281             Double_t  innerPerpSA  = tTrackSA.innermostMeasurementState().globalPosition().perp();
00282             Double_t  outerPerpSA  = tTrackSA.outermostMeasurementState().globalPosition().perp();
00283      
00284             TrajectoryStateOnSurface innerTSOS=tTrackSA.outermostMeasurementState();
00285 //      PropagationDirection propagationDir=alongMomentum;
00286             const Propagator * thePropagator;
00287       
00288             // Define which kind of reco track is used
00289             if ( (outerPerpSA-innerPerpSA) > 0 ) {
00290 
00291                 trackRefitterType = "LHCLike";
00292                 innerTSOS = tTrackSA.innermostMeasurementState();
00293                 thePropagator = thePropagatorAlo.product();
00294 //      propagationDir = alongMomentum;
00295           
00296             }else {//if ((outerPerpSA-innerPerpSA) < 0 ) {
00297         
00298                 trackRefitterType = "CosmicLike";
00299                 innerTSOS = tTrackSA.outermostMeasurementState();
00300                 thePropagator = thePropagatorOpp.product(); 
00301 //      propagationDir = oppositeToMomentum;
00302       
00303             }   
00304 
00305             RecHitVector  my4DTrack = this->doMatching(*muon, all4DSegmentsDT, all4DSegmentsCSC, &indexCollectionDT, &indexCollectionCSC, theTrackingGeometry);
00306   
00307     
00308 //cut in number of segments
00309             if(my4DTrack.size()>(min4DTrackSegmentSize-1) ){
00310 
00311 
00312 // start propagation
00313 //    TrajectoryStateOnSurface innerTSOS = track.impactPointState();
00314 //                    TrajectoryStateOnSurface innerTSOS = track.innermostMeasurementState();
00315 
00316                 //If the state is valid
00317                 if(innerTSOS.isValid()) {
00318 
00319                     //Loop over Associated segments
00320                     for(RecHitVector::iterator rechit = my4DTrack.begin(); rechit != my4DTrack.end(); ++rechit) {
00321         
00322                         const GeomDet* geomDet = theTrackingGeometry->idToDet((*rechit)->geographicalId());
00323 //Otherwise the propagator could throw an exception
00324                         const Plane* pDest = dynamic_cast<const Plane*>(&geomDet->surface());
00325                         const Cylinder* cDest = dynamic_cast<const Cylinder*>(&geomDet->surface());
00326 
00327                         if(pDest != 0 || cDest != 0) {   
00328  
00329 //                                  Propagator *updatePropagator=thePropagator->clone();
00330 //                                  updatePropagator->setPropagationDirection(propagationDir);
00331                             TrajectoryStateOnSurface destiny = thePropagator->propagate(*(innerTSOS.freeState()), geomDet->surface());
00332 
00333                             if(!destiny.isValid()|| !destiny.hasError()) continue;
00334 
00335                             const long rawId= (*rechit)->geographicalId().rawId();
00336                             int position = -1;
00337 
00338                             DetId myDet(rawId);
00339                             int det = myDet.subdetId();
00340                             int wheel=0,station=0,sector=0;
00341                             int endcap=0,ring=0,chamber=0;
00342                             bool goAhead = (det==1 && doDT) || (det==2 && doCSC);
00343 
00344                             double residualLocalX=0,residualLocalPhi=0,residualLocalY=0,residualLocalTheta=0;
00345 
00346                             // Fill generic histograms
00347                             //If it's a DT
00348                             if(det == 1 && doDT) {
00349                                 DTChamberId myChamber(rawId);
00350                                 wheel=myChamber.wheel();
00351                                 station = myChamber.station();
00352                                 sector=myChamber.sector();
00353               
00354                                 residualLocalX = (*rechit)->localPosition().x() -destiny.localPosition().x();
00355                 
00356 
00357                                 residualLocalPhi = atan2(((RecSegment *)(*rechit))->localDirection().z(), 
00358                                                          ((RecSegment*)(*rechit))->localDirection().x()) - atan2(destiny.localDirection().z(), destiny.localDirection().x());
00359                                 if(station!=4){
00360                 
00361                                     residualLocalY = (*rechit)->localPosition().y() - destiny.localPosition().y();
00362                 
00363 
00364                                     residualLocalTheta = atan2(((RecSegment *)(*rechit))->localDirection().z(), 
00365                                                                ((RecSegment*)(*rechit))->localDirection().y()) - atan2(destiny.localDirection().z(), destiny.localDirection().y());
00366 
00367                                 }
00368 
00369                             } 
00370                             else if (det==2 && doCSC){
00371                                 CSCDetId myChamber(rawId);
00372                                 endcap= myChamber.endcap();
00373                                 station = myChamber.station();
00374                                 if(endcap==2) station = -station;
00375                                 ring = myChamber.ring();
00376                                 chamber=myChamber.chamber();
00377 
00378                                 residualLocalX = (*rechit)->localPosition().x() -destiny.localPosition().x();
00379                 
00380                                 residualLocalY = (*rechit)->localPosition().y() - destiny.localPosition().y();
00381                 
00382                                 residualLocalPhi = atan2(((RecSegment *)(*rechit))->localDirection().y(), 
00383                                                          ((RecSegment*)(*rechit))->localDirection().x()) - atan2(destiny.localDirection().y(), destiny.localDirection().x());
00384 
00385                                 residualLocalTheta = atan2(((RecSegment *)(*rechit))->localDirection().y(), 
00386                                                            ((RecSegment*)(*rechit))->localDirection().z()) - atan2(destiny.localDirection().y(), destiny.localDirection().z());
00387 
00388                             }
00389                             else{
00390                                 residualLocalX=0,residualLocalPhi=0,residualLocalY=0,residualLocalTheta=0;
00391                             }
00392                         
00393                             // Fill individual chamber histograms
00394 
00395 
00396                             std::string nameOfHistoLocalX;
00397 
00398         
00399                             if(det==1 && doDT){ // DT
00400                                 std::stringstream Wheel; Wheel<<wheel;
00401                                 std::stringstream Station; Station<<station;
00402                                 std::stringstream Sector; Sector<<sector;
00403                                         
00404                                 nameOfHistoLocalX="ResidualLocalX_W"+Wheel.str()+"MB"+Station.str()+"S"+Sector.str();
00405  
00406                             }
00407                             else if(det==2 && doCSC){ //CSC
00408                                 std::stringstream Ring; Ring<<ring;
00409                                 std::stringstream Station; Station<<station;
00410                                 std::stringstream Chamber; Chamber<<chamber;
00411                                                                                
00412                                 nameOfHistoLocalX="ResidualLocalX_ME"+Station.str()+"R"+Ring.str()+"C"+Chamber.str();
00413                             }               
00414             
00415                             if (goAhead){
00416                                     
00417                                 for(unsigned int i=0 ; i<unitsLocalX.size() ; i++)
00418                                 {
00419 
00420                                     if( nameOfHistoLocalX==unitsLocalX[i]->getName()){
00421                                         position=i; break;}
00422                                 }
00423                                             
00424                                             
00425                                     if(trackRefitterType == "CosmicLike") { //problem with angle convention in reverse extrapolation
00426                                     residualLocalPhi += 3.1416;
00427                                     residualLocalTheta +=3.1416;
00428                                 }
00429                                 
00430                                 unitsLocalX.at(position)->Fill(residualLocalX);
00431                                 unitsLocalPhi.at(position)->Fill(residualLocalPhi);
00432 
00433                                 if(det==1 && station!=4) {unitsLocalY.at(position)->Fill(residualLocalY); 
00434                                     unitsLocalTheta.at(position)->Fill(residualLocalTheta);}
00435  
00436                                 else if(det==2) {unitsLocalY.at(position)->Fill(residualLocalY);
00437                                     unitsLocalTheta.at(position)->Fill(residualLocalTheta);}
00438                                 
00439            
00440                                 countPoints++;
00441                                 // if at least one point from this track is used, count this track
00442                                 if (countPoints==1) countTracks++;
00443                             }           
00444         
00445                             innerTSOS = destiny;
00446 
00447                             //delete thePropagator;
00448 
00449                         }else {
00450                             edm::LogError("MuonAlignment") <<" Error!! Exception in propagator catched" << std::endl;
00451                             continue;
00452                         }
00453 
00454                     } //loop over my4DTrack
00455                 } //TSOS was valid
00456 
00457             } // cut in at least 4 segments
00458 
00459         } //end cut in RecHitsSize>36
00460         numberOfHits=numberOfHits+countPoints;
00461 //       } //Muon is GlobalMuon
00462     } //loop over Muons
00463     numberOfTracks=numberOfTracks+countTracks;
00464 
00465 //        delete thePropagator;
00466 
00467 //      else edm::LogError("MuonAlignment")<<"Error!! Specified MuonCollection "<< theMuonTrackCollectionLabel <<" is not present in event. ProductNotFound!!"<<std::endl;
00468 }
00469 
00470 RecHitVector MuonAlignment::doMatching(const reco::Track &staTrack, edm::Handle<DTRecSegment4DCollection> &all4DSegmentsDT, edm::Handle<CSCSegmentCollection> &all4DSegmentsCSC, intDVector *indexCollectionDT, intDVector *indexCollectionCSC, edm::ESHandle<GlobalTrackingGeometry> &theTrackingGeometry) {
00471 
00472     DTRecSegment4DCollection::const_iterator segmentDT;
00473     CSCSegmentCollection::const_iterator segmentCSC;
00474   
00475     std::vector<int> positionDT;
00476     std::vector<int> positionCSC;
00477     RecHitVector my4DTrack;
00478   
00479     //Loop over the hits of the track
00480     for(unsigned int counter = 0; counter != staTrack.recHitsSize()-1; counter++) {
00481     
00482         TrackingRecHitRef myRef = staTrack.recHit(counter);
00483         const TrackingRecHit *rechit = myRef.get();
00484         const GeomDet* geomDet = theTrackingGeometry->idToDet(rechit->geographicalId());
00485     
00486         //It's a DT Hit
00487         if(geomDet->subDetector() == GeomDetEnumerators::DT) {
00488       
00489             //Take the layer associated to this hit
00490             DTLayerId myLayer(rechit->geographicalId().rawId());
00491       
00492             int NumberOfDTSegment = 0;
00493             //Loop over segments
00494             for(segmentDT = all4DSegmentsDT->begin(); segmentDT != all4DSegmentsDT->end(); ++segmentDT) {
00495         
00496                 //By default the chamber associated to this Segment is new
00497                 bool isNewChamber = true;
00498         
00499                 //Loop over segments already included in the vector of segments in the actual track
00500                 for(std::vector<int>::iterator positionIt = positionDT.begin(); positionIt != positionDT.end(); positionIt++) {
00501           
00502                     //If this segment has been used before isNewChamber = false
00503                     if(NumberOfDTSegment == *positionIt) isNewChamber = false;
00504                 }
00505         
00506                 //Loop over vectors of segments associated to previous tracks
00507                 for(std::vector<std::vector<int> >::iterator collect = indexCollectionDT->begin(); collect != indexCollectionDT->end(); ++collect) {
00508           
00509                     //Loop over segments associated to a track
00510                     for(std::vector<int>::iterator positionIt = (*collect).begin(); positionIt != (*collect).end(); positionIt++) {
00511             
00512                         //If this segment was used in a previos track then isNewChamber = false
00513                         if(NumberOfDTSegment == *positionIt) isNewChamber = false;
00514                     }
00515                 }
00516         
00517                 //If the chamber is new
00518                 if(isNewChamber) {
00519           
00520                     DTChamberId myChamber((*segmentDT).geographicalId().rawId());
00521                     //If the layer of the hit belongs to the chamber of the 4D Segment
00522                     if(myLayer.wheel() == myChamber.wheel() && myLayer.station() == myChamber.station() && myLayer.sector() == myChamber.sector()) {
00523             
00524                         //push position of the segment and tracking rechit
00525                         positionDT.push_back(NumberOfDTSegment);
00526                         my4DTrack.push_back((TrackingRecHit *) &(*segmentDT));
00527                     }
00528                 }
00529                 NumberOfDTSegment++;
00530             }
00531             //In case is a CSC
00532         } else if (geomDet->subDetector() == GeomDetEnumerators::CSC) {
00533       
00534             //Take the layer associated to this hit
00535             CSCDetId myLayer(rechit->geographicalId().rawId());
00536       
00537             int NumberOfCSCSegment = 0;
00538             //Loop over 4Dsegments
00539             for(segmentCSC = all4DSegmentsCSC->begin(); segmentCSC != all4DSegmentsCSC->end(); segmentCSC++) {
00540         
00541                 //By default the chamber associated to the segment is new
00542                 bool isNewChamber = true;
00543         
00544                 //Loop over segments in the current track
00545                 for(std::vector<int>::iterator positionIt = positionCSC.begin(); positionIt != positionCSC.end(); positionIt++) {
00546           
00547                     //If this segment has been used then newchamber = false
00548                     if(NumberOfCSCSegment == *positionIt) isNewChamber = false;
00549                 }
00550                 //Loop over vectors of segments in previous tracks
00551                 for(std::vector<std::vector<int> >::iterator collect = indexCollectionCSC->begin(); collect != indexCollectionCSC->end(); ++collect) {
00552           
00553                     //Loop over segments in a track
00554                     for(std::vector<int>::iterator positionIt = (*collect).begin(); positionIt != (*collect).end(); positionIt++) {
00555             
00556                         //If the segment was used in a previous track isNewChamber = false
00557                         if(NumberOfCSCSegment == *positionIt) isNewChamber = false;
00558                     }
00559                 }
00560                 //If the chamber is new
00561                 if(isNewChamber) {
00562           
00563                     CSCDetId myChamber((*segmentCSC).geographicalId().rawId());
00564                     //If the chambers are the same
00565                     if(myLayer.chamberId() == myChamber.chamberId()) {
00566                         //push
00567                         positionCSC.push_back(NumberOfCSCSegment);
00568                         my4DTrack.push_back((TrackingRecHit *) &(*segmentCSC));
00569                     }
00570                 }
00571                 NumberOfCSCSegment++;
00572             }
00573         }
00574     }
00575   
00576     indexCollectionDT->push_back(positionDT);
00577     indexCollectionCSC->push_back(positionCSC);
00578   
00579     if ( trackRefitterType == "CosmicLike") {
00580 
00581         std::reverse(my4DTrack.begin(),my4DTrack.end());
00582 
00583     }
00584     return my4DTrack;
00585 
00586 
00587 }
00588 
00589 
00590 
00591 void MuonAlignment::endJob(void) {
00592 
00593 
00594     LogTrace(metname)<<"[MuonAlignment] Saving the histos";
00595     bool outputMEsInRootFile = parameters.getParameter<bool>("OutputMEsInRootFile");
00596     std::string outputFileName = parameters.getParameter<std::string>("OutputFileName");
00597 
00598     edm::LogInfo("MuonAlignment")  << "Number of Tracks considered for residuals: " << numberOfTracks << std::endl << std::endl;
00599     edm::LogInfo("MuonAlignment")  << "Number of Hits considered for residuals: " << numberOfHits << std::endl << std::endl;
00600 
00601     if (doSummary){
00602         char binLabel[15];
00603 
00604         // check if ME still there (and not killed by MEtoEDM for memory saving)
00605         if( dbe )
00606           {
00607             // check existence of first histo in the list
00608             if (! dbe->get(topFolder.str()+"/DT/hLocalPositionDT")) return;
00609           }
00610         else 
00611           return;
00612         
00613 
00614         for(unsigned int i=0 ; i<unitsLocalX.size() ; i++)
00615         {
00616 
00617             if(unitsLocalX[i]->getEntries()!=0){
00618 
00619                 TString nameHistoLocalX=unitsLocalX[i]->getName();
00620 
00621                 TString nameHistoLocalPhi=unitsLocalPhi[i]->getName();
00622 
00623                 TString nameHistoLocalTheta=unitsLocalTheta[i]->getName();
00624 
00625                 TString nameHistoLocalY=unitsLocalY[i]->getName();
00626 
00627 
00628                 if (nameHistoLocalX.Contains("MB") ) // HistoLocalX DT
00629                 {
00630                     int wheel, station, sector;
00631 
00632                     sscanf(nameHistoLocalX, "ResidualLocalX_W%dMB%1dS%d",&wheel,&station,&sector);
00633 
00634                     Int_t nstation=station - 1;
00635                     Int_t nwheel=wheel+2;
00636 
00637                     Double_t Mean=unitsLocalX[i]->getMean();
00638                     Double_t Error=unitsLocalX[i]->getMeanError();
00639 
00640                     Int_t ybin=1+nwheel*8+nstation*2;
00641                     hLocalPositionDT->setBinContent(sector,ybin,fabs(Mean));
00642                     sprintf(binLabel, "MB%d/%d_X",wheel, station );
00643                     hLocalPositionDT->setBinLabel(ybin,binLabel,2);
00644                     hLocalPositionRmsDT->setBinContent(sector,ybin,Error);
00645                     hLocalPositionRmsDT->setBinLabel(ybin,binLabel,2);
00646                 
00647                     hLocalXMeanDT->Fill(fabs(Mean));
00648                     hLocalXRmsDT->Fill(Error);
00649                 }
00650 
00651                 if (nameHistoLocalX.Contains("ME")) // HistoLocalX CSC
00652                 {
00653                     int station, ring, chamber;
00654 
00655                     sscanf(nameHistoLocalX, "ResidualLocalX_ME%dR%1dC%d",&station,&ring,&chamber);
00656 
00657                     Double_t Mean=unitsLocalX[i]->getMean();
00658                     Double_t Error=unitsLocalX[i]->getMeanError();
00659 
00660                     Int_t ybin=abs(station)*2+ring;
00661                     if(abs(station)==1) ybin=ring;
00662                     if (station>0) ybin=ybin+10;
00663                     else ybin = 11 -ybin;
00664                     ybin=2*ybin-1;
00665                     hLocalPositionCSC->setBinContent(chamber,ybin,fabs(Mean));
00666                     sprintf(binLabel, "ME%d/%d_X", station,ring );
00667                     hLocalPositionCSC->setBinLabel(ybin,binLabel,2);
00668                     hLocalPositionRmsCSC->setBinContent(chamber,ybin,Error);
00669                     hLocalPositionRmsCSC->setBinLabel(ybin,binLabel,2);
00670 
00671                     hLocalXMeanCSC->Fill(fabs(Mean));
00672                     hLocalXRmsCSC->Fill(Error);
00673                 }
00674 
00675                 if (nameHistoLocalTheta.Contains("MB")) // HistoLocalTheta DT
00676                 {       
00677 
00678                     int wheel, station, sector;
00679 
00680                     sscanf(nameHistoLocalTheta, "ResidualLocalTheta_W%dMB%1dS%d",&wheel,&station,&sector);
00681 
00682                     if(station != 4){
00683                         Int_t nstation=station - 1;
00684                         Int_t nwheel=wheel+2;
00685 
00686                         Double_t Mean=unitsLocalTheta[i]->getMean();
00687                         Double_t Error=unitsLocalTheta[i]->getMeanError();
00688 
00689                         Int_t ybin=2+nwheel*8+nstation*2;
00690                         hLocalAngleDT->setBinContent(sector,ybin,fabs(Mean));
00691                         sprintf(binLabel, "MB%d/%d_#theta",wheel,station );
00692                         hLocalAngleDT->setBinLabel(ybin,binLabel,2);
00693                         hLocalAngleRmsDT->setBinContent(sector,ybin,Error);
00694                         hLocalAngleRmsDT->setBinLabel(ybin,binLabel,2);
00695         
00696                         hLocalThetaMeanDT->Fill(fabs(Mean));
00697                         hLocalThetaRmsDT->Fill(Error);
00698                     }
00699                 }
00700 
00701                 if (nameHistoLocalPhi.Contains("MB")) // HistoLocalPhi DT
00702                 {       
00703 
00704                     int wheel, station, sector;
00705 
00706                     sscanf(nameHistoLocalPhi, "ResidualLocalPhi_W%dMB%1dS%d",&wheel,&station,&sector);
00707 
00708                     Int_t nstation=station - 1;
00709                     Int_t nwheel=wheel+2;
00710 
00711                     Double_t Mean=unitsLocalPhi[i]->getMean();
00712                     Double_t Error=unitsLocalPhi[i]->getMeanError();
00713 
00714                     Int_t ybin=1+nwheel*8+nstation*2;
00715                     hLocalAngleDT->setBinContent(sector,ybin,fabs(Mean));
00716                     sprintf(binLabel, "MB%d/%d_#phi", wheel,station );
00717                     hLocalAngleDT->setBinLabel(ybin,binLabel,2);
00718                     hLocalAngleRmsDT->setBinContent(sector,ybin,Error);
00719                     hLocalAngleRmsDT->setBinLabel(ybin,binLabel,2);
00720 
00721                     hLocalPhiMeanDT->Fill(fabs(Mean));
00722                     hLocalPhiRmsDT->Fill(Error);
00723                 }
00724 
00725                 if (nameHistoLocalPhi.Contains("ME")) // HistoLocalPhi CSC
00726                 {
00727 
00728                     int station, ring, chamber;
00729 
00730                     sscanf(nameHistoLocalPhi, "ResidualLocalPhi_ME%dR%1dC%d",&station,&ring,&chamber);
00731 
00732                     Double_t Mean=unitsLocalPhi[i]->getMean();
00733                     Double_t Error=unitsLocalPhi[i]->getMeanError();
00734 
00735                     Int_t ybin=abs(station)*2+ring;
00736                     if(abs(station)==1) ybin=ring;
00737                     if (station>0) ybin=ybin+10;
00738                     else ybin = 11 -ybin;
00739                     ybin=2*ybin-1;
00740                     hLocalAngleCSC->setBinContent(chamber,ybin,fabs(Mean));
00741                     sprintf(binLabel, "ME%d/%d_#phi", station,ring );
00742                     hLocalAngleCSC->setBinLabel(ybin,binLabel,2);
00743                     hLocalAngleRmsCSC->setBinContent(chamber,ybin,Error);
00744                     hLocalAngleRmsCSC->setBinLabel(ybin,binLabel,2);
00745 
00746                     hLocalPhiMeanCSC->Fill(fabs(Mean));
00747                     hLocalPhiRmsCSC->Fill(Error);
00748                 }
00749 
00750                 if (nameHistoLocalTheta.Contains("ME")) // HistoLocalTheta CSC
00751                 {
00752 
00753                     int station, ring, chamber;
00754 
00755                     sscanf(nameHistoLocalTheta, "ResidualLocalTheta_ME%dR%1dC%d",&station,&ring,&chamber);
00756 
00757                     Double_t Mean=unitsLocalTheta[i]->getMean();
00758                     Double_t Error=unitsLocalTheta[i]->getMeanError();
00759 
00760                     Int_t ybin=abs(station)*2+ring;
00761                     if(abs(station)==1) ybin=ring;
00762                     if (station>0) ybin=ybin+10;
00763                     else ybin = 11 -ybin;
00764                     ybin=2*ybin;
00765                     hLocalAngleCSC->setBinContent(chamber,ybin,fabs(Mean));
00766                     sprintf(binLabel, "ME%d/%d_#theta", station,ring );
00767                     hLocalAngleCSC->setBinLabel(ybin,binLabel,2);
00768                     hLocalAngleRmsCSC->setBinContent(chamber,ybin,Error);
00769                     hLocalAngleRmsCSC->setBinLabel(ybin,binLabel,2);
00770 
00771                     hLocalThetaMeanCSC->Fill(fabs(Mean));
00772                     hLocalThetaRmsCSC->Fill(Error);
00773 
00774                 }
00775 
00776                 if (nameHistoLocalY.Contains("MB")) // HistoLocalY DT
00777                 {
00778 
00779                     int wheel, station, sector;
00780 
00781                     sscanf(nameHistoLocalY, "ResidualLocalY_W%dMB%1dS%d",&wheel,&station,&sector);
00782 
00783                     if(station!=4){
00784                         Int_t nstation=station - 1;
00785                         Int_t nwheel=wheel+2;
00786 
00787                         Double_t Mean=unitsLocalY[i]->getMean();
00788                         Double_t Error=unitsLocalY[i]->getMeanError();
00789 
00790                         Int_t ybin=2+nwheel*8+nstation*2;
00791                         hLocalPositionDT->setBinContent(sector,ybin,fabs(Mean));
00792                         sprintf(binLabel, "MB%d/%d_Y", wheel,station );
00793                         hLocalPositionDT->setBinLabel(ybin,binLabel,2);
00794                         hLocalPositionRmsDT->setBinContent(sector,ybin,Error);
00795                         hLocalPositionRmsDT->setBinLabel(ybin,binLabel,2);
00796 
00797                         hLocalYMeanDT->Fill(fabs(Mean));
00798                         hLocalYRmsDT->Fill(Error);
00799                     }
00800                 }
00801 
00802                 if (nameHistoLocalY.Contains("ME")) // HistoLocalY CSC
00803                 {
00804 
00805                     int station, ring, chamber;
00806 
00807                     sscanf(nameHistoLocalY, "ResidualLocalY_ME%dR%1dC%d",&station,&ring,&chamber);
00808 
00809                     Double_t Mean=unitsLocalY[i]->getMean();
00810                     Double_t Error=unitsLocalY[i]->getMeanError();
00811 
00812                     Int_t ybin=abs(station)*2+ring;
00813                     if(abs(station)==1) ybin=ring;
00814                     if (station>0) ybin=ybin+10;
00815                     else ybin = 11 -ybin;
00816                     ybin=2*ybin;
00817                     hLocalPositionCSC->setBinContent(chamber,ybin,fabs(Mean));
00818                     sprintf(binLabel, "ME%d/%d_Y", station,ring );
00819                     hLocalPositionCSC->setBinLabel(ybin,binLabel,2);
00820                     hLocalPositionRmsCSC->setBinContent(chamber,ybin,Error);
00821                     hLocalPositionRmsCSC->setBinLabel(ybin,binLabel,2);
00822 
00823                     hLocalYMeanCSC->Fill(fabs(Mean));
00824                     hLocalYRmsCSC->Fill(Error);
00825                 }
00826             } // check in # entries
00827         } // loop on vector of histos
00828     } //doSummary
00829     
00830     if(outputMEsInRootFile){
00831 //    dbe->showDirStructure();
00832         dbe->save(outputFileName);
00833     }
00834 
00835 
00836 }
00837
/data/doxygen/doxygen-1.7.3/gen/CMSSW_4_2_8/src/DQMOffline/Alignment/src/MuonAlignment.cc
