CMS 3D CMS Logo

/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_7/src/RecoLocalMuon/CSCEfficiency/src/CSCEfficiency.cc

Go to the documentation of this file.
00001 /*
00002  *  Routine to calculate CSC efficiencies 
00003  *  Comments about the program logic are denoted by //----
00004  * 
00005  *  Stoyan Stoynev, Northwestern University.
00006  */ 
00007   
00008 #include "RecoLocalMuon/CSCEfficiency/src/CSCEfficiency.h"
00009 
00010 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
00011 #include "TrackPropagation/SteppingHelixPropagator/interface/SteppingHelixPropagator.h"
00012 #include "FWCore/Common/interface/TriggerNames.h"
00013 #include "DataFormats/MuonReco/interface/Muon.h"
00014 #include "DataFormats/MuonReco/interface/MuonFwd.h"
00015 #include "DataFormats/MuonReco/interface/MuonSelectors.h"
00016 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
00017 using namespace std;
00018 using namespace edm;
00019 
00020 //---- The Analysis  (main)
00021 bool CSCEfficiency::filter(Event & event, const EventSetup& eventSetup){
00022   passTheEvent = false;
00023   DataFlow->Fill(0.);  
00024   MuonPatternRecoDumper debug;
00025  
00026   //---- increment counter
00027   nEventsAnalyzed++;
00028   // printalot debug output
00029   printalot = (nEventsAnalyzed < int(printout_NEvents)); // 
00030   int iRun   = event.id().run();
00031   int iEvent = event.id().event();
00032   if(0==fmod(double (nEventsAnalyzed) ,double(1000) )){
00033     if(printalot){
00034       printf("\n==enter==CSCEfficiency===== run %i\tevent %i\tn Analyzed %i\n",iRun,iEvent,nEventsAnalyzed);
00035     }
00036   }
00037   theService->update(eventSetup);  
00038   //---- These declarations create handles to the types of records that you want
00039   //---- to retrieve from event "e".
00040   if (printalot) printf("\tget handles for digi collections\n");
00041   
00042   //---- Pass the handle to the method "getByType", which is used to retrieve
00043   //---- one and only one instance of the type in question out of event "e". If
00044   //---- zero or more than one instance exists in the event an exception is thrown.
00045   if (printalot) printf("\tpass handles\n");
00046   edm::Handle<CSCALCTDigiCollection> alcts;
00047   edm::Handle<CSCCLCTDigiCollection> clcts;
00048   edm::Handle<CSCCorrelatedLCTDigiCollection> correlatedlcts;
00049   edm::Handle<CSCWireDigiCollection> wires;
00050   edm::Handle<CSCStripDigiCollection> strips;
00051   edm::Handle<CSCRecHit2DCollection> rechits; 
00052   edm::Handle<CSCSegmentCollection> segments;
00053   //edm::Handle<reco::TrackCollection> saMuons;
00054   edm::Handle<edm::View<reco::Track> > trackCollectionH;
00055   edm::Handle<edm::PSimHitContainer> simhits;
00056 
00057   if(useDigis){
00058     event.getByLabel(alctDigiTag_, alcts);
00059     event.getByLabel(clctDigiTag_, clcts);
00060     event.getByLabel(corrlctDigiTag_, correlatedlcts);
00061 
00062     event.getByLabel( stripDigiTag_, strips);
00063     event.getByLabel( wireDigiTag_,  wires);
00064   }
00065   if(!isData){
00066     event.getByLabel(simHitTag, simhits);
00067   }
00068   event.getByLabel(rechitDigiTag_,rechits); 
00069   event.getByLabel(segmentDigiTag_, segments);
00070   //event.getByLabel(saMuonTag,saMuons);
00071   event.getByLabel(tracksTag,trackCollectionH);
00072   const edm::View<reco::Track>  trackCollection = *(trackCollectionH.product());
00073 
00074   //---- Get the CSC Geometry :
00075   if (printalot) printf("\tget the CSC geometry.\n");
00076   edm::ESHandle<CSCGeometry> cscGeom;
00077   eventSetup.get<MuonGeometryRecord>().get(cscGeom);
00078 
00079   // use theTrackingGeometry instead of cscGeom?
00080   ESHandle<GlobalTrackingGeometry> theTrackingGeometry;
00081   eventSetup.get<GlobalTrackingGeometryRecord>().get(theTrackingGeometry);
00082 
00083   bool triggerPassed = true;
00084   if(useTrigger){
00085     // access the trigger information
00086     // trigger names can be find in HLTrigger/Configuration/python/HLT_2E30_cff.py (or?)
00087    // get hold of TriggerResults
00088     edm::Handle<edm::TriggerResults> hltR;
00089     event.getByLabel(hlTriggerResults_,hltR);
00090     const edm::TriggerNames & triggerNames = event.triggerNames(*hltR);
00091     triggerPassed = applyTrigger(hltR, triggerNames);
00092   }
00093   if(!triggerPassed){
00094     return triggerPassed;
00095   }
00096   DataFlow->Fill(1.);   
00097   GlobalPoint gpZero(0.,0.,0.);
00098   if(theService->magneticField()->inTesla(gpZero).mag2()<0.1){ 
00099     magField = false;
00100   }
00101   else{
00102     magField = true;
00103   }
00104 
00105   //---- store info from digis
00106   fillDigiInfo(alcts, clcts, correlatedlcts, wires, strips, simhits, rechits, segments, cscGeom);
00107   //
00108   Handle<reco::MuonCollection> muons;
00109   edm::InputTag muonTag_("muons");
00110   event.getByLabel(muonTag_,muons);
00111 
00112   edm::Handle<reco::BeamSpot> beamSpotHandle;
00113   event.getByLabel("offlineBeamSpot", beamSpotHandle);
00114   reco::BeamSpot vertexBeamSpot = *beamSpotHandle;
00115   //
00116   std::vector <reco::MuonCollection::const_iterator> goodMuons_it;
00117   unsigned int nPositiveZ = 0;
00118   unsigned int nNegativeZ = 0;
00119   float muonOuterZPosition = -99999.;
00120   if(isIPdata){
00121     if (printalot)std::cout<<" muons.size() = "<<muons->size() <<std::endl;
00122     for ( reco::MuonCollection::const_iterator muon = muons->begin(); muon != muons->end(); ++muon ) {
00123       DataFlow->Fill(31.);      
00124       if (printalot) {
00125         std::cout<<"  iMuon = "<<muon-muons->begin()<<" charge = "<<muon->charge()<<" p = "<<muon->p()<<" pt = "<<muon->pt()<<
00126           " eta = "<<muon->eta()<<" phi = "<<muon->phi()<<
00127           " matches = "<<
00128           muon->matches().size()<<" matched Seg = "<<muon->numberOfMatches(reco::Muon::SegmentAndTrackArbitration)<<" GLB/TR/STA = "<<
00129           muon->isGlobalMuon()<<"/"<<muon->isTrackerMuon()<<"/"<<muon->isStandAloneMuon()<<std::endl;
00130       }
00131       if(!(muon->isTrackerMuon() && muon->isGlobalMuon())){
00132         continue;
00133       }
00134       DataFlow->Fill(32.);
00135       double relISO = ( muon->isolationR03().sumPt +
00136                         muon->isolationR03().emEt +
00137                         muon->isolationR03().hadEt)/muon->track()->pt();
00138       if (printalot) {
00139         std::cout<<" relISO = "<<relISO<<" emVetoEt = "<<muon->isolationR03().emVetoEt<<" caloComp = "<<
00140           muon::caloCompatibility(*(muon))<<" dxy = "<<fabs(muon->track()->dxy(vertexBeamSpot.position()))<<std::endl;
00141       }
00142       if(
00143          //relISO>0.1 || muon::caloCompatibility(*(muon))<.90 ||
00144          fabs(muon->track()->dxy(vertexBeamSpot.position()))>0.2 || muon->pt()<6.){
00145         continue;
00146       }
00147       DataFlow->Fill(33.);
00148       if(muon->track()->hitPattern().numberOfValidPixelHits()<1 ||
00149          muon->track()->hitPattern().numberOfValidTrackerHits()<11 ||
00150          muon->combinedMuon()->hitPattern().numberOfValidMuonHits()<1 ||
00151          muon->combinedMuon()->normalizedChi2()>10. ||
00152          muon->numberOfMatches()<2){
00153         continue;
00154       }
00155       DataFlow->Fill(34.);
00156       float zOuter = muon->combinedMuon()->outerPosition().z();
00157       float rhoOuter = muon->combinedMuon()->outerPosition().rho();
00158       bool passDepth = true;
00159       // barrel region 
00160       //if ( fabs(zOuter) < 660. && rhoOuter > 400. && rhoOuter < 480.){
00161       if ( fabs(zOuter) < 660. && rhoOuter > 400. && rhoOuter < 540.){ 
00162         passDepth = false;
00163       }
00164       // endcap region
00165       //else if( fabs(zOuter) > 550. && fabs(zOuter) < 650. && rhoOuter < 300.){
00166       else if( fabs(zOuter) > 550. && fabs(zOuter) < 650. && rhoOuter < 300.){
00167         passDepth = false;
00168       }
00169       // overlap region
00170       //else if ( fabs(zOuter) > 680. && fabs(zOuter) < 730. && rhoOuter < 480.){
00171       else if ( fabs(zOuter) > 680. && fabs(zOuter) < 880. && rhoOuter < 540.){
00172         passDepth = false;
00173       }
00174       if(!passDepth){
00175         continue;
00176       }
00177       DataFlow->Fill(35.);
00178       goodMuons_it.push_back(muon);
00179       if(muon->track()->momentum().z()>0.){
00180         ++nPositiveZ;
00181       }
00182       if(muon->track()->momentum().z()<0.){
00183         ++nNegativeZ;
00184       }
00185     }
00186   }
00187 
00188   //
00189 
00190 
00191   if (printalot) std::cout<<"Start track loop over "<<trackCollection.size()<<" tracks"<<std::endl;
00192   for(edm::View<reco::Track>::size_type i=0; i<trackCollection.size(); ++i) {
00193     DataFlow->Fill(2.);
00194     edm::RefToBase<reco::Track> track(trackCollectionH, i);
00195     //std::cout<<" iTR = "<<i<<" eta = "<<track->eta()<<" phi = "<<track->phi()<<std::cout<<" pt = "<<track->pt()<<std::endl;
00196     if(isIPdata){
00197       if (printalot){
00198         std::cout<<" nNegativeZ = "<<nNegativeZ<<" nPositiveZ = "<<nPositiveZ<<std::endl;
00199       }
00200       if(nNegativeZ>1 || nPositiveZ>1){
00201         break;
00202       }
00203       bool trackOK = false;
00204       if (printalot){
00205         std::cout<<" goodMuons_it.size() = "<<goodMuons_it.size()<<std::endl;
00206       }
00207       for(size_t iM=0;iM<goodMuons_it.size();++iM){
00208         //std::cout<<" iM = "<<iM<<" eta = "<<goodMuons_it[iM]->track()->eta()<<
00209         //" phi = "<<goodMuons_it[iM]->track()->phi()<<
00210         //" pt = "<<goodMuons_it[iM]->track()->pt()<<std::endl;
00211         float deltaR = pow(track->phi()-goodMuons_it[iM]->track()->phi(),2) +
00212           pow(track->eta()-goodMuons_it[iM]->track()->eta(),2);
00213         deltaR = sqrt(deltaR);
00214         if (printalot){
00215           std::cout<<" TR mu match to a tr: deltaR = "<<deltaR<<" dPt = "<<
00216             track->pt()-goodMuons_it[iM]->track()->pt()<<std::endl;
00217         }
00218         if(deltaR>0.01 || fabs(track->pt()-goodMuons_it[iM]->track()->pt())>0.1 ){
00219           continue;
00220         }
00221         else{
00222           trackOK = true;
00223           if (printalot){
00224             std::cout<<" trackOK "<<std::endl;
00225           }
00226           muonOuterZPosition = goodMuons_it[iM]->combinedMuon()->outerPosition().z();
00227           break;
00228           //++nChosenTracks;
00229         }
00230       }
00231       if(!trackOK){
00232         if (printalot){
00233           std::cout<<" failed: trackOK "<<std::endl;
00234         }
00235         continue;
00236       }
00237     }
00238     else{
00239       //---- Do we need a better "clean track" definition?
00240       if(trackCollection.size()>2){
00241         break;
00242       }
00243       DataFlow->Fill(3.);
00244       if(!i && 2==trackCollection.size()){
00245         edm::View<reco::Track>::size_type tType = 1;
00246         edm::RefToBase<reco::Track> trackTwo(trackCollectionH, tType);
00247         if(track->outerPosition().z()*trackTwo->outerPosition().z()>0){// in one and the same "endcap"
00248           break;
00249         }
00250       }
00251     }
00252     DataFlow->Fill(4.);
00253     if (printalot){
00254       std::cout<<"i track = "<<i<<" P = "<<track->p()<<" chi2/ndf = "<<track->normalizedChi2()<<" nSeg = "<<segments->size()<<std::endl;
00255       std::cout<<"quality undef/loose/tight/high/confirmed/goodIt/size "<<
00256         track->quality(reco::Track::undefQuality)<<"/"<<
00257         track->quality(reco::Track::loose)<<"/"<<
00258         track->quality(reco::Track::tight)<<"/"<<
00259         track->quality(reco::Track::highPurity)<<"/"<<
00260         track->quality(reco::Track::confirmed)<<"/"<<
00261         track->quality(reco::Track::goodIterative)<<"/"<<
00262         track->quality(reco::Track::qualitySize)<<
00263         std::endl;
00264       std::cout<<" pt = "<< track->pt()<<" +-"<<track->ptError()<<" q/pt = "<<track->qoverp()<<" +- "<<track->qoverpError()<<std::endl;
00265       //std::cout<<" const Pmin = "<<minTrackMomentum<<" pMax = "<<maxTrackMomentum<<" maxNormChi2 = "<<maxNormChi2<<std::endl;
00266       std::cout<<" track inner position = "<<track->innerPosition()<<" outer position = "<<track->outerPosition()<<std::endl;
00267       std::cout<<"track eta (outer) = "<<track->outerPosition().eta()<<" phi (outer) = "<<
00268         track->outerPosition().phi()<<std::endl;
00269       if(fabs(track->innerPosition().z())>500.){      
00270         DetId innerDetId(track->innerDetId());
00271         std::cout<<" dump inner state MUON detid  = "<<debug.dumpMuonId(innerDetId)<<std::endl;
00272       }
00273       if(fabs(track->outerPosition().z())>500.){
00274         DetId outerDetId(track->outerDetId());
00275         std::cout<<" dump outer state MUON detid  = "<<debug.dumpMuonId(outerDetId)<<std::endl;
00276       }
00277 
00278       std::cout<<" nHits = "<<track->found()<<std::endl;
00279       /*
00280       trackingRecHit_iterator rhbegin = track->recHitsBegin();
00281       trackingRecHit_iterator rhend = track->recHitsEnd();
00282       int iRH = 0;
00283       for(trackingRecHit_iterator recHit = rhbegin; recHit != rhend; ++recHit){
00284         const GeomDet* geomDet = theTrackingGeometry->idToDet((*recHit)->geographicalId());
00285         std::cout<<"hit "<<iRH<<" loc pos = " <<(*recHit)->localPosition()<<
00286           " glob pos = " <<geomDet->toGlobal((*recHit)->localPosition())<<std::endl;
00287         ++iRH;
00288       }
00289       */
00290     }
00291     float dpT_ov_pT = 0.;
00292     if(fabs(track->pt())>0.001){
00293       dpT_ov_pT =  track->ptError()/ track->pt();
00294     }
00295     //---- These define a "good" track
00296     if(track->normalizedChi2()>maxNormChi2){// quality
00297       break;
00298     }
00299     DataFlow->Fill(5.);
00300     if(track->found()<minTrackHits){// enough data points
00301       break;
00302     }
00303     DataFlow->Fill(6.);
00304     if(!segments->size()){// better have something in the CSC 
00305       break;
00306     }
00307     DataFlow->Fill(7.);
00308     if(magField && (track->p()<minP || track->p()>maxP)){// proper energy range 
00309       break;
00310     }
00311     DataFlow->Fill(8.);
00312     if(magField && (dpT_ov_pT >0.5) ){// not too crazy uncertainty
00313       break;
00314     }
00315     DataFlow->Fill(9.);
00316 
00317     passTheEvent = true; 
00318     if (printalot) std::cout<<"good Track"<<std::endl;
00319     CLHEP::Hep3Vector r3T_inner(track->innerPosition().x(),track->innerPosition().y(),track->innerPosition().z());
00320     CLHEP::Hep3Vector r3T(track->outerPosition().x(),track->outerPosition().y(),track->outerPosition().z());
00321     chooseDirection(r3T_inner, r3T);// for non-IP
00322 
00323     CLHEP::Hep3Vector p3T(track->outerMomentum().x(),track->outerMomentum().y(),track->outerMomentum().z());
00324     CLHEP::Hep3Vector p3_propagated, r3_propagated;
00325     AlgebraicSymMatrix66 cov_propagated, covT;
00326     covT *= 1e-20;
00327     cov_propagated *= 1e-20;
00328     int charge = track->charge(); 
00329     FreeTrajectoryState ftsStart = getFromCLHEP(p3T, r3T, charge, covT, &*(theService->magneticField()));
00330     if (printalot){
00331       std::cout<<" p = "<<track->p()<<" norm chi2 = "<<track->normalizedChi2()<<std::endl;
00332       std::cout<<" dump the very first FTS  = "<<debug.dumpFTS(ftsStart)<<std::endl;
00333     }
00334     TrajectoryStateOnSurface tSOSDest;
00335     int endcap = 0;
00336     //---- which endcap to look at
00337     if(track->outerPosition().z()>0){
00338       endcap = 1;
00339     }
00340     else{
00341       endcap = 2;
00342     }
00343     int chamber = 1;
00344     //---- a "refference" CSCDetId for each ring
00345     std::vector< CSCDetId > refME;
00346     for(int iS=1;iS<5;++iS){
00347       for(int iR=1;iR<4;++iR){
00348         if(1!=iS && iR>2){
00349           continue;
00350         }
00351         else if(4==iS && iR>1){
00352           continue;
00353         } 
00354         refME.push_back( CSCDetId(endcap, iS, iR, chamber));
00355       }
00356     }
00357     //---- loop over the "refference" CSCDetIds
00358     for(size_t iSt = 0; iSt<refME.size();++iSt){
00359       if (printalot){
00360         std::cout<<"loop iStatation = "<<iSt<<std::endl;
00361         std::cout<<"refME[iSt]: st = "<<refME[iSt].station()<<" rg = "<<refME[iSt].ring()<<std::endl;
00362       }
00363       std::map <std::string, bool> chamberTypes;
00364       chamberTypes["ME11"] = false;
00365       chamberTypes["ME12"] = false;
00366       chamberTypes["ME13"] = false;
00367       chamberTypes["ME21"] = false;
00368       chamberTypes["ME22"] = false;
00369       chamberTypes["ME31"] = false;
00370       chamberTypes["ME32"] = false;
00371       chamberTypes["ME41"] = false;
00372       const CSCChamber* cscChamber_base = cscGeom->chamber(refME[iSt].chamberId());
00373       DetId detId = cscChamber_base->geographicalId();
00374       if (printalot){
00375         std::cout<<" base iStation : eta = "<<cscGeom->idToDet(detId)->surface().position().eta()<<" phi = "<<
00376           cscGeom->idToDet(detId)->surface().position().phi() << " y = " <<cscGeom->idToDet(detId)->surface().position().y()<<std::endl;
00377         std::cout<<" dump base iStation detid  = "<<debug.dumpMuonId(detId)<<std::endl;
00378         std::cout<<" dump FTS start  = "<<debug.dumpFTS(ftsStart)<<std::endl;
00379       }
00380       //---- propagate to this ME
00381       tSOSDest = propagate(ftsStart, cscGeom->idToDet(detId)->surface());
00382       if(tSOSDest.isValid()){
00383         ftsStart = *tSOSDest.freeState();
00384         if (printalot) std::cout<<"  dump FTS end   = "<<debug.dumpFTS(ftsStart)<<std::endl;
00385         getFromFTS(ftsStart, p3_propagated, r3_propagated, charge, cov_propagated);
00386         float feta =  fabs(r3_propagated.eta());
00387         float phi =  r3_propagated.phi();
00388         //---- which rings are (possibly) penetrated
00389         ringCandidates(refME[iSt].station(), feta, chamberTypes);
00390 
00391         map<std::string,bool>::iterator iter;   
00392         int iterations = 0;
00393         //---- loop over ring candidates 
00394         for( iter = chamberTypes.begin(); iter != chamberTypes.end(); iter++ ) {
00395           ++iterations;
00396           //---- is this ME a machinig candidate station 
00397           if(iter->second && (iterations-1)==int(iSt)){
00398             if (printalot){
00399               std::cout<<" Chamber type "<< iter->first<<" is a candidate..."<<std::endl;
00400               std::cout<<" station() = "<< refME[iSt].station()<<" ring() = "<<refME[iSt].ring()<<" iSt = "<<iSt<<std::endl;
00401             }
00402             std::vector <int> coupleOfChambers;
00403             //---- which chamber (and its closes neighbor) is penetrated by the track - candidates
00404             chamberCandidates(refME[iSt].station(), refME[iSt].ring(), phi, coupleOfChambers);
00405             //---- loop over the two chamber candidates
00406             for(size_t iCh =0;iCh<coupleOfChambers.size();++iCh){
00407               DataFlow->Fill(11.);  
00408               if (printalot) std::cout<<" Check chamber N = "<<coupleOfChambers.at(iCh)<<std::endl;;
00409               if((!getAbsoluteEfficiency) && (true == emptyChambers
00410                                               [refME[iSt].endcap()-1]
00411                                               [refME[iSt].station()-1]
00412                                               [refME[iSt].ring()-1]
00413                                               [coupleOfChambers.at(iCh)-FirstCh])){
00414                 continue;
00415               }
00416               CSCDetId theCSCId(refME[iSt].endcap(), refME[iSt].station(), refME[iSt].ring(), coupleOfChambers.at(iCh));
00417               const CSCChamber* cscChamber = cscGeom->chamber(theCSCId.chamberId());
00418               const BoundPlane bpCh = cscGeom->idToDet(cscChamber->geographicalId())->surface();
00419               float zFTS = ftsStart.position().z();
00420               float dz = fabs(bpCh.position().z() - zFTS);
00421               float zDistInner = track->innerPosition().z() - bpCh.position().z();
00422               float zDistOuter = track->outerPosition().z() - bpCh.position().z();
00423               //---- only detectors between the inner and outer points of the track are considered for non IP-data
00424               if(printalot){
00425                 std::cout<<" zIn = "<<track->innerPosition().z()<<" zOut = "<<track->outerPosition().z()<<" zSurf = "<<bpCh.position().z()<<std::endl;
00426               }
00427               if(!isIPdata && (zDistInner*zDistOuter>0. || fabs(zDistInner)<15. || fabs(zDistOuter)<15.)){ // for non IP-data
00428                 if(printalot){
00429                   std::cout<<" Not an intermediate (as defined) point... Skip."<<std::endl;
00430                 }
00431                 continue;
00432               }
00433               if(isIPdata && fabs(track->eta())<1.8){
00434                 if(fabs(muonOuterZPosition) - fabs(bpCh.position().z())<0 || 
00435                    fabs(muonOuterZPosition-bpCh.position().z())<15.){
00436                   continue;
00437                 }
00438               }
00439               DataFlow->Fill(13.);
00440               //---- propagate to the chamber (from this ME) if it is a different surface (odd/even chambers)
00441               if(dz>0.1){// i.e. non-zero (float 0 check is bad) 
00442                 //if(fabs(zChanmber - zFTS ) > 0.1){
00443                 tSOSDest = propagate(ftsStart, cscGeom->idToDet(cscChamber->geographicalId())->surface());                   
00444                 if(tSOSDest.isValid()){
00445                   ftsStart = *tSOSDest.freeState();
00446                 }
00447                 else{
00448                   if(printalot) std::cout<<"TSOS not valid! Break."<<std::endl;
00449                   break;
00450                 }
00451               }
00452               else{
00453                 if(printalot) std::cout<<" info: dz<0.1"<<std::endl;
00454               }
00455               DataFlow->Fill(15.);  
00456               FreeTrajectoryState ftsInit = ftsStart;
00457               bool inDeadZone = false;
00458               //---- loop over the 6 layers
00459               for(int iLayer = 0;iLayer<6;++iLayer){
00460                 bool extrapolationPassed = true;
00461                 if (printalot){
00462                   std::cout<<" iLayer = "<<iLayer<<"   dump FTS init  = "<<debug.dumpFTS(ftsInit)<<std::endl;
00463                   std::cout<<" dump detid  = "<<debug.dumpMuonId(cscChamber->geographicalId())<<std::endl;
00464                   std::cout<<"Surface to propagate to:  pos = "<<cscChamber->layer(iLayer+1)->surface().position()<<" eta = "
00465                            <<cscChamber->layer(iLayer+1)->surface().position().eta()<<" phi = "
00466                            <<cscChamber->layer(iLayer+1)->surface().position().phi()<<std::endl;
00467                 }
00468                 //---- propagate to this layer
00469                 tSOSDest = propagate(ftsInit, cscChamber->layer(iLayer+1)->surface());
00470                 if(tSOSDest.isValid()){
00471                   ftsInit = *tSOSDest.freeState();
00472                   if (printalot) std::cout<<" Propagation between layers successful:  dump FTS end  = "<<debug.dumpFTS(ftsInit)<<std::endl;
00473                   getFromFTS(ftsInit, p3_propagated, r3_propagated, charge, cov_propagated);
00474                 }
00475                 else{
00476                   if (printalot) std::cout<<"Propagation between layers not successful - notValid TSOS"<<std::endl;
00477                   extrapolationPassed = false;
00478                   inDeadZone = true;
00479                 }
00480                 //}
00481                 //---- Extrapolation passed? For each layer?
00482                 if(extrapolationPassed){
00483                   GlobalPoint theExtrapolationPoint(r3_propagated.x(),r3_propagated.y(),r3_propagated.z());
00484                   LocalPoint theLocalPoint = cscChamber->layer(iLayer+1)->toLocal(theExtrapolationPoint);
00485                   //std::cout<<" Candidate chamber: extrapolated LocalPoint = "<<theLocalPoint<<std::endl;
00486                   inDeadZone = ( inDeadZone ||
00487                                  !inSensitiveLocalRegion(theLocalPoint.x(), theLocalPoint.y(), 
00488                                                          refME[iSt].station(), refME[iSt].ring()));
00489                   if (printalot){
00490                     std::cout<<" Candidate chamber: extrapolated LocalPoint = "<<theLocalPoint<<"inDeadZone = "<<inDeadZone<<std::endl;
00491                   }
00492                   //---- break if in dead zone for any layer ("clean" tracks)
00493                   if(inDeadZone){
00494                     break;
00495                   }
00496                 }
00497                 else{
00498                   break;
00499                 }
00500               }
00501               DataFlow->Fill(17.);
00502               //---- Is a track in a sensitive area for each layer?  
00503               if(!inDeadZone){//---- for any layer
00504                 DataFlow->Fill(19.);  
00505                 if (printalot) std::cout<<"Do efficiencies..."<<std::endl;
00506                 //---- Do efficiencies
00507                 // angle cuts applied (if configured)
00508                 bool angle_flag = true; angle_flag = efficienciesPerChamber(theCSCId, cscChamber, ftsStart);  
00509                 if(useDigis && angle_flag){
00510                   stripWire_Efficiencies(theCSCId, ftsStart);
00511                 }
00512                 if(angle_flag){
00513                   recHitSegment_Efficiencies(theCSCId, cscChamber, ftsStart);
00514                   if(!isData){
00515                     recSimHitEfficiency(theCSCId, ftsStart);
00516                   }
00517                 }
00518               }
00519               else{
00520                 if(printalot) std::cout<<" Not in active area for all layers"<<std::endl;
00521               }
00522             }
00523             if(tSOSDest.isValid()){
00524               ftsStart = *tSOSDest.freeState();
00525             }
00526           }
00527         }
00528       }
00529       else{
00530         if (printalot) std::cout<<" TSOS not valid..."<<std::endl;
00531       }
00532     }
00533   }
00534   //---- End
00535   if (printalot) printf("==exit===CSCEfficiency===== run %i\tevent %i\n\n",iRun,iEvent);
00536   return passTheEvent;
00537 }
00538 
00539 //
00540 bool CSCEfficiency::inSensitiveLocalRegion(double xLocal, double yLocal, int station, int ring){
00541   //---- Good region means sensitive area of a chamber. "Local" stands for the local system 
00542   bool pass = false;
00543   std::vector <double> chamberBounds(3);// the sensitive area
00544   float y_center = 99999.;
00545   //---- hardcoded... not good
00546   if(station>1 && station<5){
00547     if(2==ring){
00548       chamberBounds[0] = 66.46/2; // (+-)x1 shorter
00549       chamberBounds[1] = 127.15/2; // (+-)x2 longer 
00550       chamberBounds[2] = 323.06/2;
00551       y_center = -0.95;
00552     }
00553     else{
00554       if(2==station){
00555         chamberBounds[0] = 54.00/2; // (+-)x1 shorter
00556         chamberBounds[1] = 125.71/2; // (+-)x2 longer 
00557         chamberBounds[2] = 189.66/2;
00558         y_center = -0.955;
00559       }
00560       else if(3==station){
00561         chamberBounds[0] = 61.40/2; // (+-)x1 shorter
00562         chamberBounds[1] = 125.71/2; // (+-)x2 longer 
00563         chamberBounds[2] = 169.70/2;
00564         y_center = -0.97;
00565       }
00566       else if(4==station){
00567         chamberBounds[0] = 69.01/2; // (+-)x1 shorter
00568         chamberBounds[1] = 125.65/2; // (+-)x2 longer 
00569         chamberBounds[2] = 149.42/2;
00570         y_center = -0.94;
00571       }
00572     }
00573   }
00574   else if(1==station){
00575     if(3==ring){
00576       chamberBounds[0] = 63.40/2; // (+-)x1 shorter
00577       chamberBounds[1] = 92.10/2; // (+-)x2 longer 
00578       chamberBounds[2] = 164.16/2;
00579       y_center = -1.075;
00580     }
00581     else if(2==ring){
00582       chamberBounds[0] = 51.00/2; // (+-)x1 shorter
00583       chamberBounds[1] = 83.74/2; // (+-)x2 longer 
00584       chamberBounds[2] = 174.49/2;
00585       y_center = -0.96;
00586     }
00587     else{// to be investigated
00588       chamberBounds[0] = 30./2;//40./2; // (+-)x1 shorter
00589       chamberBounds[1] = 60./2;//100./2; // (+-)x2 longer 
00590       chamberBounds[2] = 160./2;//142./2;
00591       y_center = 0.;
00592     }
00593   }
00594   double yUp = chamberBounds[2] + y_center;
00595   double yDown = - chamberBounds[2] + y_center;
00596   double xBound1Shifted = chamberBounds[0]-distanceFromDeadZone;//
00597   double xBound2Shifted = chamberBounds[1]-distanceFromDeadZone;//
00598   double lineSlope = (yUp - yDown)/(xBound2Shifted-xBound1Shifted);
00599   double lineConst = yUp - lineSlope*xBound2Shifted;
00600   double yBoundary =  lineSlope*abs(xLocal) + lineConst;
00601   pass = checkLocal(yLocal, yBoundary, station, ring);
00602   return pass;
00603 }
00604 
00605 bool CSCEfficiency::checkLocal(double yLocal, double yBoundary, int station, int ring){
00606 //---- check if it is in a good local region (sensitive area - geometrical and HV boundaries excluded) 
00607   bool pass = false;
00608   std::vector <float> deadZoneCenter(6);
00609   const float deadZoneHalf = 0.32*7/2;// wire spacing * (wires missing + 1)/2
00610   float cutZone = deadZoneHalf + distanceFromDeadZone;//cm
00611   //---- hardcoded... not good
00612   if(station>1 && station<5){
00613     if(2==ring){
00614       deadZoneCenter[0]= -162.48 ;
00615       deadZoneCenter[1] = -81.8744;
00616       deadZoneCenter[2] = -21.18165;
00617       deadZoneCenter[3] = 39.51105;
00618       deadZoneCenter[4] = 100.2939;
00619       deadZoneCenter[5] = 160.58;
00620       
00621       if(yLocal >yBoundary &&
00622          ((yLocal> deadZoneCenter[0] + cutZone && yLocal< deadZoneCenter[1] - cutZone) ||
00623           (yLocal> deadZoneCenter[1] + cutZone && yLocal< deadZoneCenter[2] - cutZone) ||
00624           (yLocal> deadZoneCenter[2] + cutZone && yLocal< deadZoneCenter[3] - cutZone) ||
00625           (yLocal> deadZoneCenter[3] + cutZone && yLocal< deadZoneCenter[4] - cutZone) ||
00626           (yLocal> deadZoneCenter[4] + cutZone && yLocal< deadZoneCenter[5] - cutZone))){
00627         pass = true;
00628       }
00629     }
00630     else if(1==ring){
00631       if(2==station){
00632         deadZoneCenter[0]= -95.94 ;
00633         deadZoneCenter[1] = -27.47;
00634         deadZoneCenter[2] = 33.67;
00635         deadZoneCenter[3] = 93.72;
00636       }
00637       else if(3==station){
00638         deadZoneCenter[0]= -85.97 ;
00639         deadZoneCenter[1] = -36.21;
00640         deadZoneCenter[2] = 23.68;
00641         deadZoneCenter[3] = 84.04;
00642       }
00643       else if(4==station){
00644         deadZoneCenter[0]= -75.82;
00645         deadZoneCenter[1] = -26.14;
00646         deadZoneCenter[2] = 23.85;
00647         deadZoneCenter[3] = 73.91;
00648       }
00649       if(yLocal >yBoundary &&
00650          ((yLocal> deadZoneCenter[0] + cutZone && yLocal< deadZoneCenter[1] - cutZone) ||
00651           (yLocal> deadZoneCenter[1] + cutZone && yLocal< deadZoneCenter[2] - cutZone) ||
00652           (yLocal> deadZoneCenter[2] + cutZone && yLocal< deadZoneCenter[3] - cutZone))){
00653         pass = true;
00654       }
00655     }
00656   }
00657   else if(1==station){
00658     if(3==ring){
00659       deadZoneCenter[0]= -83.155 ;
00660       deadZoneCenter[1] = -22.7401;
00661       deadZoneCenter[2] = 27.86665;
00662       deadZoneCenter[3] = 81.005;
00663       if(yLocal > yBoundary &&
00664          ((yLocal> deadZoneCenter[0] + cutZone && yLocal< deadZoneCenter[1] - cutZone) ||
00665           (yLocal> deadZoneCenter[1] + cutZone && yLocal< deadZoneCenter[2] - cutZone) ||
00666           (yLocal> deadZoneCenter[2] + cutZone && yLocal< deadZoneCenter[3] - cutZone))){
00667         pass = true;
00668       }
00669     }
00670     else if(2==ring){
00671       deadZoneCenter[0]= -86.285 ;
00672       deadZoneCenter[1] = -32.88305;
00673       deadZoneCenter[2] = 32.867423;
00674       deadZoneCenter[3] = 88.205;
00675       if(yLocal > (yBoundary) &&
00676          ((yLocal> deadZoneCenter[0] + cutZone && yLocal< deadZoneCenter[1] - cutZone) ||
00677           (yLocal> deadZoneCenter[1] + cutZone && yLocal< deadZoneCenter[2] - cutZone) ||
00678           (yLocal> deadZoneCenter[2] + cutZone && yLocal< deadZoneCenter[3] - cutZone))){
00679         pass = true;
00680       }
00681     }
00682     else{
00683       deadZoneCenter[0]= -81.0;
00684       deadZoneCenter[1] = 81.0;
00685       if(yLocal > (yBoundary) &&
00686          ((yLocal> deadZoneCenter[0] + cutZone && yLocal< deadZoneCenter[1] - cutZone) )){
00687         pass = true;
00688       }
00689     }
00690   }
00691   return pass;
00692 }
00693 
00694 void CSCEfficiency::fillDigiInfo(edm::Handle<CSCALCTDigiCollection> &alcts, 
00695                                  edm::Handle<CSCCLCTDigiCollection> &clcts, 
00696                                  edm::Handle<CSCCorrelatedLCTDigiCollection> &correlatedlcts,
00697                                  edm::Handle<CSCWireDigiCollection> &wires,
00698                                  edm::Handle<CSCStripDigiCollection> &strips,
00699                                  edm::Handle<edm::PSimHitContainer> &simhits,
00700                                  edm::Handle<CSCRecHit2DCollection> &rechits,
00701                                  edm::Handle<CSCSegmentCollection> &segments,
00702                                  edm::ESHandle<CSCGeometry> &cscGeom){
00703   for(int iE=0;iE<2;iE++){
00704     for(int iS=0;iS<4;iS++){
00705       for(int iR=0;iR<4;iR++){
00706         for(int iC=0;iC<NumCh;iC++){
00707           allSegments[iE][iS][iR][iC].clear(); 
00708           allCLCT[iE][iS][iR][iC] = allALCT[iE][iS][iR][iC] = allCorrLCT[iE][iS][iR][iC] = false;
00709           for(int iL=0;iL<6;iL++){
00710             allStrips[iE][iS][iR][iC][iL].clear(); 
00711             allWG[iE][iS][iR][iC][iL].clear();
00712             allRechits[iE][iS][iR][iC][iL].clear(); 
00713             allSimhits[iE][iS][iR][iC][iL].clear(); 
00714           }
00715         }
00716       }
00717     }
00718   }
00719   //
00720   if(useDigis){
00721     fillLCT_info(alcts, clcts, correlatedlcts);
00722     fillWG_info(wires, cscGeom);
00723     fillStrips_info(strips);
00724   }
00725   fillRechitsSegments_info(rechits, segments, cscGeom);
00726   if(!isData){
00727     fillSimhit_info(simhits);
00728   }
00729 }
00730 
00731 
00732 void CSCEfficiency::fillLCT_info(edm::Handle<CSCALCTDigiCollection> &alcts, 
00733                                  edm::Handle<CSCCLCTDigiCollection> &clcts, 
00734                                  edm::Handle<CSCCorrelatedLCTDigiCollection> &correlatedlcts ){
00735   //---- ALCTDigis
00736   int nSize = 0;
00737   for (CSCALCTDigiCollection::DigiRangeIterator j=alcts->begin(); j!=alcts->end(); j++) {
00738     ++nSize;
00739     const CSCDetId& id = (*j).first;
00740     const CSCALCTDigiCollection::Range& range =(*j).second;
00741     for (CSCALCTDigiCollection::const_iterator digiIt =
00742            range.first; digiIt!=range.second;
00743          ++digiIt){
00744       // Valid digi in the chamber (or in neighbouring chamber) 
00745       if((*digiIt).isValid()){
00746         allALCT[id.endcap()-1][id.station()-1][id.ring()-1][id.chamber()-FirstCh] = true;
00747       }
00748     }// for digis in layer
00749   }// end of for (j=...
00750   ALCTPerEvent->Fill(nSize);
00751   //---- CLCTDigis
00752   nSize = 0;
00753   for (CSCCLCTDigiCollection::DigiRangeIterator j=clcts->begin(); j!=clcts->end(); j++) {
00754     ++nSize;
00755     const CSCDetId& id = (*j).first;
00756     std::vector<CSCCLCTDigi>::const_iterator digiIt = (*j).second.first;
00757     std::vector<CSCCLCTDigi>::const_iterator last = (*j).second.second;
00758     for( ; digiIt != last; ++digiIt) {
00759       // Valid digi in the chamber (or in neighbouring chamber) 
00760       if((*digiIt).isValid()){
00761         allCLCT[id.endcap()-1][id.station()-1][id.ring()-1][id.chamber()-FirstCh] = true;
00762       }
00763     }
00764   }
00765   CLCTPerEvent->Fill(nSize);
00766   //---- CorrLCTDigis
00767   for (CSCCorrelatedLCTDigiCollection::DigiRangeIterator j=correlatedlcts->begin(); j!=correlatedlcts->end(); j++) {
00768     const CSCDetId& id = (*j).first;
00769     std::vector<CSCCorrelatedLCTDigi>::const_iterator digiIt = (*j).second.first;
00770     std::vector<CSCCorrelatedLCTDigi>::const_iterator last = (*j).second.second;
00771     for( ; digiIt != last; ++digiIt) {
00772       // Valid digi in the chamber (or in neighbouring chamber) 
00773       if((*digiIt).isValid()){
00774         allCorrLCT[id.endcap()-1][id.station()-1][id.ring()-1][id.chamber()-FirstCh] = true;
00775       }
00776     }
00777   }
00778 }
00779 //
00780 void CSCEfficiency::fillWG_info(edm::Handle<CSCWireDigiCollection> &wires, edm::ESHandle<CSCGeometry> &cscGeom){
00781   //---- WIRE GROUPS
00782   for (CSCWireDigiCollection::DigiRangeIterator j=wires->begin(); j!=wires->end(); j++) {
00783     CSCDetId id = (CSCDetId)(*j).first;
00784     const CSCLayer *layer_p = cscGeom->layer (id);
00785     const CSCLayerGeometry *layerGeom = layer_p->geometry ();
00786     //
00787     std::vector<CSCWireDigi>::const_iterator digiItr = (*j).second.first;
00788     std::vector<CSCWireDigi>::const_iterator last = (*j).second.second;
00789     //
00790     for( ; digiItr != last; ++digiItr) {
00791       std::pair < int, float > WG_pos(digiItr->getWireGroup(), layerGeom->yOfWireGroup(digiItr->getWireGroup())); 
00792       std::pair <std::pair < int, float >, int >  LayerSignal(WG_pos, digiItr->getTimeBin()); 
00793       
00794       //---- AllWG contains basic information about WG (WG number and Y-position, time bin)
00795       allWG[id.endcap()-1][id.station()-1][id.ring()-1][id.chamber()-FirstCh]
00796         [id.layer()-1].push_back(LayerSignal);
00797       if(printalot){
00798         //std::cout<<" WG check : "<<std::endl;
00799         //printf("\t\tendcap/station/ring/chamber/layer: %i/%i/%i/%i/%i\n",id.endcap(),id.station(),id.ring(),id.chamber(),id.layer());
00800         //std::cout<<" WG size = "<<allWG[id.endcap()-1][id.station()-1][id.ring()-1][id.chamber()-FirstCh]
00801         //[id.layer()-1].size()<<std::endl;
00802       }
00803     }
00804   }
00805 }
00806 void CSCEfficiency::fillStrips_info(edm::Handle<CSCStripDigiCollection> &strips){
00807   //---- STRIPS
00808   for (CSCStripDigiCollection::DigiRangeIterator j=strips->begin(); j!=strips->end(); j++) {
00809     CSCDetId id = (CSCDetId)(*j).first;
00810     int largestADCValue = -1;
00811     std::vector<CSCStripDigi>::const_iterator digiItr = (*j).second.first;
00812     std::vector<CSCStripDigi>::const_iterator last = (*j).second.second;
00813     for( ; digiItr != last; ++digiItr) {
00814       int maxADC=largestADCValue;
00815       int myStrip = digiItr->getStrip();
00816       std::vector<int> myADCVals = digiItr->getADCCounts();
00817       float thisPedestal = 0.5*(float)(myADCVals[0]+myADCVals[1]);
00818       float threshold = 13.3 ;
00819       float diff = 0.;
00820       float peakADC  = -1000.;
00821       for (unsigned int iCount = 0; iCount < myADCVals.size(); iCount++) {
00822         diff = (float)myADCVals[iCount]-thisPedestal;
00823         if (diff > threshold) { 
00824           if (myADCVals[iCount] > largestADCValue) {
00825             largestADCValue = myADCVals[iCount];
00826           }
00827         }
00828         if (diff > threshold && diff > peakADC) {
00829           peakADC  = diff;
00830         }
00831       }
00832       if(largestADCValue>maxADC){// FIX IT!!!
00833         maxADC = largestADCValue;
00834         std::pair <int, float> LayerSignal (myStrip, peakADC);
00835         
00836         //---- AllStrips contains basic information about strips 
00837         //---- (strip number and peak signal for most significant strip in the layer) 
00838         allStrips[id.endcap()-1][id.station()-1][id.ring()-1][id.chamber()-1][id.layer()-1].clear();
00839         allStrips[id.endcap()-1][id.station()-1][id.ring()-1][id.chamber()-1][id.layer()-1].push_back(LayerSignal);
00840       }
00841     }
00842   }
00843 }
00844 void CSCEfficiency::fillSimhit_info(edm::Handle<edm::PSimHitContainer> &simhits){
00845   //---- SIMHITS
00846   edm::PSimHitContainer::const_iterator dSHsimIter;
00847   for (dSHsimIter = simhits->begin(); dSHsimIter != simhits->end(); dSHsimIter++){
00848     // Get DetID for this simHit:
00849     CSCDetId sId = (CSCDetId)(*dSHsimIter).detUnitId();
00850     std::pair <LocalPoint, int> simHitPos((*dSHsimIter).localPosition(), (*dSHsimIter).particleType());  
00851     allSimhits[sId.endcap()-1][sId.station()-1][sId.ring()-1][sId.chamber()-FirstCh][sId.layer()-1].push_back(simHitPos);
00852   }
00853 }
00854 //
00855 void CSCEfficiency::fillRechitsSegments_info(edm::Handle<CSCRecHit2DCollection> &rechits,
00856                                              edm::Handle<CSCSegmentCollection> &segments,
00857                                              edm::ESHandle<CSCGeometry> &cscGeom
00858                                              ){
00859   //---- RECHITS AND SEGMENTS
00860   //---- Loop over rechits 
00861   if (printalot){ 
00862     //printf("\tGet the recHits collection.\t ");
00863     printf("  The size of the rechit collection is %i\n",int(rechits->size()));
00864     //printf("\t...start loop over rechits...\n");
00865   }
00866   recHitsPerEvent->Fill(rechits->size());
00867   //---- Build iterator for rechits and loop :
00868   CSCRecHit2DCollection::const_iterator recIt;
00869   for (recIt = rechits->begin(); recIt != rechits->end(); recIt++) {
00870     //---- Find chamber with rechits in CSC 
00871     CSCDetId id = (CSCDetId)(*recIt).cscDetId();
00872     if (printalot){
00873       const CSCLayer* csclayer = cscGeom->layer( id);
00874       LocalPoint rhitlocal = (*recIt).localPosition();  
00875       LocalError rerrlocal = (*recIt).localPositionError();  
00876       GlobalPoint rhitglobal= csclayer->toGlobal(rhitlocal);
00877       printf("\t\tendcap/station/ring/chamber/layer: %i/%i/%i/%i/%i\n",id.endcap(),id.station(),id.ring(),id.chamber(),id.layer());
00878       printf("\t\tx,y,z: %f, %f, %f\texx,eey,exy: %f, %f, %f\tglobal x,y,z: %f, %f, %f \n",
00879              rhitlocal.x(), rhitlocal.y(), rhitlocal.z(), rerrlocal.xx(), rerrlocal.yy(), rerrlocal.xy(),
00880              rhitglobal.x(), rhitglobal.y(), rhitglobal.z());
00881     }
00882     std::pair <LocalPoint, bool> recHitPos((*recIt).localPosition(), false);  
00883     allRechits[id.endcap()-1][id.station()-1][id.ring()-1][id.chamber()-FirstCh][id.layer()-1].push_back(recHitPos);
00884   }
00885   //---- "Empty" chambers
00886   for(int iE=0;iE<2;iE++){
00887     for(int iS=0;iS<4;iS++){
00888       for(int iR=0;iR<4;iR++){
00889         for(int iC=0;iC<NumCh;iC++){
00890           int numLayers = 0;
00891           for(int iL=0;iL<6;iL++){
00892             if(allRechits[iE][iS][iR][iC][iL].size()){
00893               ++numLayers; 
00894             }
00895           }
00896           if(numLayers>1){
00897             emptyChambers[iE][iS][iR][iC] = false;
00898           }
00899           else{
00900             emptyChambers[iE][iS][iR][iC] = true;
00901           }
00902         }
00903       }
00904     }
00905   }
00906 
00907   //
00908   if (printalot){
00909     printf("  The size of the segment collection is %i\n", int(segments->size()));
00910     //printf("\t...start loop over segments...\n");
00911   }
00912   segmentsPerEvent->Fill(segments->size());
00913   for(CSCSegmentCollection::const_iterator it = segments->begin(); it != segments->end(); it++) {
00914     CSCDetId id  = (CSCDetId)(*it).cscDetId();
00915     StHist[id.endcap()-1][id.station()-1].segmentChi2_ndf->Fill((*it).chi2()/(*it).degreesOfFreedom());
00916     StHist[id.endcap()-1][id.station()-1].hitsInSegment->Fill((*it).nRecHits());
00917     if (printalot){ 
00918       printf("\tendcap/station/ring/chamber: %i %i %i %i\n",
00919              id.endcap(),id.station(),id.ring(),id.chamber());
00920       std::cout<<"\tposition(loc) = "<<(*it).localPosition()<<" error(loc) = "<<(*it).localPositionError()<<std::endl;
00921       std::cout<<"\t chi2/ndf = "<<(*it).chi2()/(*it).degreesOfFreedom()<<" nhits = "<<(*it).nRecHits() <<std::endl;
00922 
00923     }
00924     allSegments[id.endcap()-1][id.station()-1][id.ring()-1][id.chamber()-FirstCh].push_back
00925       (make_pair((*it).localPosition(), (*it).localDirection()));
00926 
00927 
00928     //---- try to get the CSC recHits that contribute to this segment.
00929     //if (printalot) printf("\tGet the recHits for this segment.\t");
00930     std::vector<CSCRecHit2D> theseRecHits = (*it).specificRecHits();
00931     int nRH = (*it).nRecHits();
00932     if (printalot){
00933       printf("\tGet the recHits for this segment.\t");
00934       printf("    nRH = %i\n",nRH);
00935     }
00936     //---- Find which of the rechits in the chamber is in the segment
00937     int layerRH = 0;
00938     for ( vector<CSCRecHit2D>::const_iterator iRH = theseRecHits.begin(); iRH != theseRecHits.end(); iRH++) {
00939       ++layerRH;
00940       CSCDetId idRH = (CSCDetId)(*iRH).cscDetId();
00941       if(printalot){ 
00942         printf("\t%i RH\tendcap/station/ring/chamber/layer: %i/%i/%i/%i/%i\n",
00943                layerRH,idRH.endcap(),idRH.station(),idRH.ring(),idRH.chamber(),idRH.layer());
00944       }
00945       for(size_t jRH = 0; 
00946           jRH<allRechits[idRH.endcap()-1][idRH.station()-1][idRH.ring()-1][idRH.chamber()-FirstCh][idRH.layer()-1].size();
00947           ++jRH){
00948         allRechits[idRH.endcap()-1][idRH.station()-1][idRH.ring()-1][idRH.chamber()-FirstCh][idRH.layer()-1][jRH].first;
00949         float xDiff = iRH->localPosition().x() - 
00950           allRechits[idRH.endcap()-1][idRH.station()-1][idRH.ring()-1][idRH.chamber()-FirstCh][idRH.layer()-1][jRH].first.x();
00951         float yDiff = iRH->localPosition().y() - 
00952           allRechits[idRH.endcap()-1][idRH.station()-1][idRH.ring()-1][idRH.chamber()-FirstCh][idRH.layer()-1][jRH].first.y();
00953         if(fabs(xDiff)<0.0001 && fabs(yDiff)<0.0001){
00954           std::pair <LocalPoint, bool> 
00955             recHitPos(allRechits[idRH.endcap()-1][idRH.station()-1][idRH.ring()-1][idRH.chamber()-FirstCh][idRH.layer()-1][jRH].first, true);
00956           allRechits[idRH.endcap()-1][idRH.station()-1][idRH.ring()-1][idRH.chamber()-FirstCh][idRH.layer()-1][jRH] = recHitPos;
00957           if(printalot){ 
00958             std::cout<<" number of the rechit (from zero) in the segment = "<< jRH<<std::endl;
00959           }
00960         }
00961       }
00962     }
00963   }
00964 }
00965 //
00966 
00967 void CSCEfficiency::ringCandidates(int station, float feta, std::map <std::string, bool> & chamberTypes){
00968   // yeah, hardcoded again...
00969   switch (station){
00970   case 1:
00971     if(feta>0.85 && feta<1.18){//ME13
00972       chamberTypes["ME13"] = true;
00973     }      
00974     if(feta>1.18 && feta<1.7){//ME12
00975       chamberTypes["ME12"] = true;
00976     }
00977     if(feta>1.5 && feta<2.45){//ME11
00978       chamberTypes["ME11"] = true;
00979     }
00980     break;
00981   case 2:
00982     if(feta>0.95 && feta<1.6){//ME22
00983       chamberTypes["ME22"] = true;
00984       
00985     }  
00986     if(feta>1.55 && feta<2.45){//ME21
00987       chamberTypes["ME21"] = true;
00988     }
00989     break;
00990   case 3:
00991     if(feta>1.08 && feta<1.72){//ME32
00992       chamberTypes["ME32"] = true;
00993       
00994     }  
00995     if(feta>1.69 && feta<2.45){//ME31
00996       chamberTypes["ME31"] = true;
00997     }
00998     break;
00999   case 4:
01000     if(feta>1.78 && feta<2.45){//ME41
01001       chamberTypes["ME41"] = true;
01002     }
01003     break;
01004   default:
01005     break;
01006   }
01007 }
01008 //
01009 void CSCEfficiency::chamberCandidates(int station, int ring, float phi, std::vector <int> &coupleOfChambers){
01010   coupleOfChambers.clear();
01011   // -pi< phi<+pi
01012   float phi_zero = 0.;// check! the phi at the "edge" of Ch 1
01013   float phi_const = 2.*M_PI/36.;
01014   int last_chamber = 36;
01015   int first_chamber = 1;
01016   if(1 != station && 1==ring){ // 18 chambers in the ring
01017     phi_const*=2;
01018     last_chamber /= 2;
01019   }
01020   if(phi<0.){
01021     if (printalot) std::cout<<" info: negative phi = "<<phi<<std::endl;
01022     phi += 2*M_PI;
01023   }
01024   float chamber_float = (phi - phi_zero)/phi_const;
01025   int chamber_int = int(chamber_float);
01026   if (chamber_float - float(chamber_int) -0.5 <0.){
01027     if(0!=chamber_int ){
01028       coupleOfChambers.push_back(chamber_int);
01029     }
01030     else{
01031       coupleOfChambers.push_back(last_chamber);
01032     }
01033     coupleOfChambers.push_back(chamber_int+1); 
01034     
01035   }
01036   else{
01037     coupleOfChambers.push_back(chamber_int+1);
01038     if(last_chamber!=chamber_int+1){
01039       coupleOfChambers.push_back(chamber_int+2);
01040     } 
01041     else{
01042       coupleOfChambers.push_back(first_chamber);
01043     }
01044   }
01045   if (printalot) std::cout<<" phi = "<<phi<<" phi_zero = "<<phi_zero<<" phi_const = "<<phi_const<<
01046     " candidate chambers: first ch = "<<coupleOfChambers[0]<<" second ch = "<<coupleOfChambers[1]<<std::endl;
01047 }
01048 
01049 //
01050 bool CSCEfficiency::efficienciesPerChamber(CSCDetId & id, const CSCChamber* cscChamber, FreeTrajectoryState &ftsChamber){
01051   int ec, st, rg, ch, secondRing;
01052   returnTypes(id, ec, st, rg, ch, secondRing);
01053 
01054   LocalVector localDir = cscChamber->toLocal(ftsChamber.momentum());
01055   if(printalot){
01056     std::cout<<" global dir = "<<ftsChamber.momentum()<<std::endl;
01057     std::cout<<" local dir = "<<localDir<<std::endl;
01058     std::cout<<" local theta = "<<localDir.theta()<<std::endl;
01059   }
01060   float dxdz = localDir.x()/localDir.z();
01061   float dydz = localDir.y()/localDir.z();
01062   if(2==st || 3==st){
01063     if(printalot){
01064       std::cout<<"st 3 or 4 ... flip dy/dz"<<std::endl;
01065     }
01066     dydz = - dydz;
01067   }
01068   if(printalot){
01069     std::cout<<"dy/dz = "<<dydz<<std::endl;
01070   }
01071   // Apply angle cut
01072   bool out = true;
01073   if(applyIPangleCuts){
01074     if(dydz>local_DY_DZ_Max || dydz<local_DY_DZ_Min || fabs(dxdz)>local_DX_DZ_Max){ 
01075       out = false;
01076     }
01077   }
01078 
01079   // Segments
01080   bool firstCondition = allSegments[ec][st][rg][ch].size() ? true : false;
01081   bool secondCondition = false;
01082   //---- ME1 is special as usual - ME1a and ME1b are actually one chamber
01083   if(secondRing>-1){
01084     secondCondition = allSegments[ec][st][secondRing][ch].size() ? true : false;
01085   }
01086   if(firstCondition || secondCondition){
01087     if(out){
01088       ChHist[ec][st][rg][ch].digiAppearanceCount->Fill(1);
01089     }  
01090   }
01091   else{
01092     if(out){
01093       ChHist[ec][st][rg][ch].digiAppearanceCount->Fill(0);
01094     } 
01095   }
01096 
01097   if(useDigis){
01098     // ALCTs
01099     firstCondition = allALCT[ec][st][rg][ch];
01100     secondCondition = false;
01101     if(secondRing>-1){
01102       secondCondition = allALCT[ec][st][secondRing][ch];
01103     }
01104     if(firstCondition || secondCondition){
01105       if(out){
01106         ChHist[ec][st][rg][ch].digiAppearanceCount->Fill(3);
01107       }
01108       // always apply partial angle cuts for this kind of histos 
01109       if(fabs(dxdz)<local_DX_DZ_Max){
01110         StHist[ec][st].EfficientALCT_momTheta->Fill(ftsChamber.momentum().theta());
01111         ChHist[ec][st][rg][ch].EfficientALCT_dydz->Fill(dydz);
01112       }
01113     }
01114     else{
01115       if(out){
01116         ChHist[ec][st][rg][ch].digiAppearanceCount->Fill(2);
01117       }
01118       if(fabs(dxdz)<local_DX_DZ_Max){ 
01119         StHist[ec][st].InefficientALCT_momTheta->Fill(ftsChamber.momentum().theta());
01120         ChHist[ec][st][rg][ch].InefficientALCT_dydz->Fill(dydz);
01121       }
01122       if(printalot){
01123         std::cout<<" missing ALCT (dy/dz = "<<dydz<<")";
01124         printf("\t\tendcap/station/ring/chamber: %i/%i/%i/%i\n",ec+1,st+1,rg+1,ch+1);
01125       }
01126     }
01127     
01128     // CLCTs
01129     firstCondition = allCLCT[ec][st][rg][ch];
01130     secondCondition = false;
01131     if(secondRing>-1){
01132       secondCondition = allCLCT[ec][st][secondRing][ch];
01133     }
01134     if(firstCondition || secondCondition){
01135       if(out){
01136         ChHist[ec][st][rg][ch].digiAppearanceCount->Fill(5);
01137       }
01138       if(dydz<local_DY_DZ_Max && dydz>local_DY_DZ_Min){
01139         StHist[ec][st].EfficientCLCT_momPhi->Fill(ftsChamber.momentum().phi() );// - phi chamber...
01140         ChHist[ec][st][rg][ch].EfficientCLCT_dxdz->Fill(dxdz);
01141       }
01142     }
01143     else{
01144       if(out){
01145         ChHist[ec][st][rg][ch].digiAppearanceCount->Fill(4);
01146       }
01147       if(dydz<local_DY_DZ_Max && dydz>local_DY_DZ_Min){ 
01148         StHist[ec][st].InefficientCLCT_momPhi->Fill(ftsChamber.momentum().phi());// - phi chamber...
01149         ChHist[ec][st][rg][ch].InefficientCLCT_dxdz->Fill(dxdz);
01150       }
01151       if(printalot){
01152         std::cout<<" missing CLCT  (dx/dz = "<<dxdz<<")";
01153         printf("\t\tendcap/station/ring/chamber: %i/%i/%i/%i\n",ec+1,st+1,rg+1,ch+1);
01154       }
01155     }
01156     if(out){
01157       // CorrLCTs
01158       firstCondition = allCorrLCT[ec][st][rg][ch];
01159       secondCondition = false;
01160       if(secondRing>-1){
01161         secondCondition = allCorrLCT[ec][st][secondRing][ch];
01162       }
01163       if(firstCondition || secondCondition){
01164         ChHist[ec][st][rg][ch].digiAppearanceCount->Fill(7);
01165       }
01166       else{
01167         ChHist[ec][st][rg][ch].digiAppearanceCount->Fill(6);
01168       }
01169     }
01170   }
01171   return out;
01172 }
01173 
01174 //
01175 bool CSCEfficiency::stripWire_Efficiencies(CSCDetId & id, FreeTrajectoryState &ftsChamber){
01176   int ec, st, rg, ch, secondRing;
01177   returnTypes(id, ec, st, rg, ch, secondRing);
01178 
01179   bool firstCondition, secondCondition;
01180   int missingLayers_s = 0;
01181   int missingLayers_wg = 0;
01182   for(int iLayer=0;iLayer<6;iLayer++){
01183  //----Strips
01184     if(printalot){ 
01185       printf("\t%i swEff: \tendcap/station/ring/chamber/layer: %i/%i/%i/%i/%i\n",
01186              iLayer + 1,id.endcap(),id.station(),id.ring(),id.chamber(),iLayer+1);
01187       std::cout<<" size S = "<<allStrips[id.endcap()-1][id.station()-1][id.ring()-1][id.chamber()-FirstCh][iLayer].size()<<
01188         "size W = "<<allWG[id.endcap()-1][id.station()-1][id.ring()-1][id.chamber()-FirstCh][iLayer].size()<<std::endl;
01189 
01190     }
01191     firstCondition = allStrips[ec][st][rg][ch][iLayer].size() ? true : false;
01192     //allSegments[ec][st][rg][ch].size() ? true : false;
01193     secondCondition = false;
01194     if(secondRing>-1){
01195       secondCondition = allStrips[ec][st][secondRing][ch][iLayer].size() ? true : false;
01196     }
01197     if(firstCondition || secondCondition){
01198       ChHist[ec][st][rg][ch].EfficientStrips->Fill(iLayer+1);
01199     }
01200     else{
01201       if(printalot){ 
01202         std::cout<<"missing strips ";
01203         printf("\t\tendcap/station/ring/chamber/layer: %i/%i/%i/%i/%i\n",id.endcap(),id.station(),id.ring(),id.chamber(),iLayer+1);
01204       }
01205     }
01206     // Wires
01207     firstCondition = allWG[ec][st][rg][ch][iLayer].size() ? true : false;
01208     secondCondition = false;
01209     if(secondRing>-1){
01210       secondCondition = allWG[ec][st][secondRing][ch][iLayer].size() ? true : false;
01211     }
01212     if(firstCondition || secondCondition){
01213       ChHist[ec][st][rg][ch].EfficientWireGroups->Fill(iLayer+1);
01214     }
01215     else{
01216       if(printalot){ 
01217         std::cout<<"missing wires ";
01218         printf("\t\tendcap/station/ring/chamber/layer: %i/%i/%i/%i/%i\n",id.endcap(),id.station(),id.ring(),id.chamber(),iLayer+1);
01219       }
01220     }
01221   }
01222   // Normalization
01223   if(6!=missingLayers_s){
01224     ChHist[ec][st][rg][ch].EfficientStrips->Fill(8);
01225   }
01226   if(6!=missingLayers_wg){
01227     ChHist[ec][st][rg][ch].EfficientWireGroups->Fill(8);
01228   }
01229   ChHist[ec][st][rg][ch].EfficientStrips->Fill(9);
01230   ChHist[ec][st][rg][ch].EfficientWireGroups->Fill(9);
01231 //
01232   ChHist[ec][st][rg][ch].StripWiresCorrelations->Fill(1);
01233   if(missingLayers_s!=missingLayers_wg){
01234     ChHist[ec][st][rg][ch].StripWiresCorrelations->Fill(2);
01235     if(6==missingLayers_wg){
01236       ChHist[ec][st][rg][ch].StripWiresCorrelations->Fill(3);
01237       ChHist[ec][st][rg][ch].NoWires_momTheta->Fill(ftsChamber.momentum().theta());
01238     }
01239     if(6==missingLayers_s){
01240       ChHist[ec][st][rg][ch].StripWiresCorrelations->Fill(4);
01241       ChHist[ec][st][rg][ch].NoStrips_momPhi->Fill(ftsChamber.momentum().theta());
01242     }
01243   }
01244   else if(6==missingLayers_s){
01245     ChHist[ec][st][rg][ch].StripWiresCorrelations->Fill(5);
01246   }
01247 
01248   return true;  
01249 }
01250 //
01251 bool CSCEfficiency::recSimHitEfficiency(CSCDetId & id, FreeTrajectoryState &ftsChamber){
01252   int ec, st, rg, ch, secondRing;
01253   returnTypes(id, ec, st, rg, ch, secondRing);
01254   bool firstCondition, secondCondition;
01255   for(int iLayer=0; iLayer<6;iLayer++){
01256     firstCondition = allSimhits[ec][st][rg][ch][iLayer].size() ? true : false;
01257     secondCondition = false;
01258     int thisRing = rg;
01259     if(secondRing>-1){
01260       secondCondition = allSimhits[ec][st][secondRing][ch][iLayer].size() ? true : false;
01261       if(secondCondition){
01262         thisRing = secondRing;
01263       }
01264     }
01265     if(firstCondition || secondCondition){
01266       for(size_t iSH=0;
01267           iSH<allSimhits[ec][st][thisRing][ch][iLayer].size();
01268           ++iSH){
01269         if(13 ==
01270            fabs(allSimhits[ec][st][thisRing][ch][iLayer][iSH].second)){
01271           ChHist[ec][st][rg][ch].SimSimhits->Fill(iLayer+1);
01272           if(allRechits[ec][st][thisRing][ch][iLayer].size()){
01273             ChHist[ec][st][rg][ch].SimRechits->Fill(iLayer+1);
01274           }
01275           break;
01276         }
01277       }
01278       //---- Next is not too usefull... 
01279       /*
01280       for(unsigned int iSimHits=0;
01281           iSimHits<allSimhits[id.endcap()-1][id.station()-1][id.ring()-1][id.chamber()-FirstCh][iLayer].size();
01282           iSimHits++){
01283         ChHist[ec][st][rg][id.chamber()-FirstCh].SimSimhits_each->Fill(iLayer+1);
01284       }
01285       for(unsigned int iRecHits=0;
01286           iRecHits<allRechits[id.endcap()-1][id.station()-1][id.ring()-1][id.chamber()-FirstCh][iLayer].size();
01287           iRecHits++){
01288         ChHist[ec][st][rg][id.chamber()-FirstCh].SimRechits_each->Fill(iLayer+1);
01289       }
01290       */
01291       //
01292     }   
01293   }
01294   return true;
01295 }
01296 
01297 //
01298 bool CSCEfficiency::recHitSegment_Efficiencies(CSCDetId & id, const CSCChamber* cscChamber, FreeTrajectoryState &ftsChamber){
01299   int ec, st, rg, ch, secondRing;
01300   returnTypes(id, ec, st, rg, ch, secondRing);
01301   bool firstCondition, secondCondition;
01302 
01303   std::vector <bool> missingLayers_rh(6);
01304   std::vector <int> usedInSegment(6);
01305   // Rechits
01306   if(printalot) std::cout<<"RecHits eff"<<std::endl;
01307   for(int iLayer=0;iLayer<6;++iLayer){
01308     firstCondition = allRechits[ec][st][rg][ch][iLayer].size() ? true : false;
01309     secondCondition = false;
01310     int thisRing = rg;
01311     if(secondRing>-1){
01312       secondCondition = allRechits[ec][st][secondRing][ch][iLayer].size() ? true : false;
01313       if(secondCondition){
01314         thisRing = secondRing;
01315       }
01316     }
01317     if(firstCondition || secondCondition){
01318       ChHist[ec][st][rg][ch].EfficientRechits_good->Fill(iLayer+1);
01319       for(size_t iR=0;
01320           iR<allRechits[ec][st][thisRing][ch][iLayer].size();
01321           ++iR){
01322         if(allRechits[ec][st][thisRing][ch][iLayer][iR].second){
01323           usedInSegment[iLayer] = 1;
01324           break;
01325         }
01326         else{
01327           usedInSegment[iLayer] = -1;
01328         }
01329       }
01330     }
01331     else{
01332       missingLayers_rh[iLayer] = true;
01333       if(printalot){
01334         std::cout<<"missing rechits ";
01335         printf("\t\tendcap/station/ring/chamber/layer: %i/%i/%i/%i/%i\n",id.endcap(),id.station(),id.ring(),id.chamber(),iLayer+1);
01336       }
01337     }
01338   }
01339   GlobalVector globalDir;
01340   GlobalPoint globalPos;
01341  // Segments
01342   firstCondition = allSegments[ec][st][rg][ch].size() ? true : false;
01343   secondCondition = false;
01344   int secondSize = 0;
01345   int thisRing = rg;
01346   if(secondRing>-1){
01347     secondCondition = allSegments[ec][st][secondRing][ch].size() ? true : false;
01348     secondSize = allSegments[ec][st][secondRing][ch].size();
01349     if(secondCondition){
01350       thisRing = secondRing;
01351     }
01352   }
01353   if(firstCondition || secondCondition){
01354     if (printalot) std::cout<<"segments - start ec = "<<ec<<" st = "<<st<<" rg = "<<rg<<" ch = "<<ch<<std::endl;
01355     StHist[ec][st].EfficientSegments_XY->Fill(ftsChamber.position().x(),ftsChamber.position().y());
01356     if(1==allSegments[ec][st][rg][ch].size() + secondSize){
01357       globalDir = cscChamber->toGlobal(allSegments[ec][st][thisRing][ch][0].second);
01358       globalPos = cscChamber->toGlobal(allSegments[ec][st][thisRing][ch][0].first);
01359       StHist[ec][st].EfficientSegments_eta->Fill(fabs(ftsChamber.position().eta()));
01360       double residual = sqrt(pow(ftsChamber.position().x() - globalPos.x(),2)+
01361                              pow(ftsChamber.position().y() - globalPos.y(),2)+
01362                              pow(ftsChamber.position().z() - globalPos.z(),2));
01363       if (printalot) std::cout<<" fts.position() = "<<ftsChamber.position()<<" segPos = "<<globalPos<<" res = "<<residual<< std::endl;
01364       StHist[ec][st].ResidualSegments->Fill(residual);
01365     }
01366     for(int iLayer=0;iLayer<6;++iLayer){
01367       if(printalot) std::cout<<" iLayer = "<<iLayer<<" usedInSegment = "<<usedInSegment[iLayer]<<std::endl;
01368       if(0!=usedInSegment[iLayer]){
01369         if(-1==usedInSegment[iLayer]){
01370           ChHist[ec][st][rg][ch].InefficientSingleHits->Fill(iLayer+1);
01371         }
01372         ChHist[ec][st][rg][ch].AllSingleHits->Fill(iLayer+1);
01373       }
01374       firstCondition = allRechits[ec][st][rg][ch][iLayer].size() ? true : false;
01375       secondCondition = false;
01376       if(secondRing>-1){
01377         secondCondition = allRechits[ec][st][secondRing][ch][iLayer].size() ? true : false;
01378       }
01379       float stripAngle = 99999.;
01380       std::vector<float> posXY(2);
01381       bool oneSegment = false;
01382       if(1==allSegments[ec][st][rg][ch].size() + secondSize){
01383         oneSegment = true;
01384         const BoundPlane bp = cscChamber->layer(iLayer+1)->surface();
01385         linearExtrapolation(globalPos,globalDir, bp.position().z(), posXY);
01386         GlobalPoint gp_extrapol( posXY.at(0), posXY.at(1),bp.position().z());
01387         const LocalPoint lp_extrapol = cscChamber->layer(iLayer+1)->toLocal(gp_extrapol);
01388         posXY.at(0) = lp_extrapol.x();
01389         posXY.at(1) = lp_extrapol.y();
01390         int nearestStrip = cscChamber->layer(iLayer+1)->geometry()->nearestStrip(lp_extrapol);
01391         stripAngle = cscChamber->layer(iLayer+1)->geometry()->stripAngle(nearestStrip) - M_PI/2. ;
01392       }
01393       if(firstCondition || secondCondition){
01394         ChHist[ec][st][rg][ch].EfficientRechits_inSegment->Fill(iLayer+1);
01395         if(oneSegment){
01396           ChHist[ec][st][rg][ch].Y_EfficientRecHits_inSegment[iLayer]->Fill(posXY.at(1));
01397           ChHist[ec][st][rg][ch].Phi_EfficientRecHits_inSegment[iLayer]->Fill(stripAngle);
01398         }
01399       }
01400       else{
01401         if(oneSegment){
01402           ChHist[ec][st][rg][ch].Y_InefficientRecHits_inSegment[iLayer]->Fill(posXY.at(1));
01403           ChHist[ec][st][rg][ch].Phi_InefficientRecHits_inSegment[iLayer]->Fill(stripAngle);
01404         }
01405       }
01406     }
01407   }
01408   else{
01409     StHist[ec][st].InefficientSegments_XY->Fill(ftsChamber.position().x(),ftsChamber.position().y());
01410     if(printalot){
01411       std::cout<<"missing segment "<<std::endl;
01412       printf("\t\tendcap/station/ring/chamber: %i/%i/%i/%i\n",id.endcap(),id.station(),id.ring(),id.chamber());
01413       std::cout<<" fts.position() = "<<ftsChamber.position()<<std::endl;
01414     }
01415   }
01416   // Normalization
01417   ChHist[ec][st][rg][ch].EfficientRechits_good->Fill(8);
01418   if(allSegments[ec][st][rg][ch].size()+secondSize<2){
01419     StHist[ec][st].AllSegments_eta->Fill(fabs(ftsChamber.position().eta()));
01420   }
01421   ChHist[ec][st][rg][id.chamber()-FirstCh].EfficientRechits_inSegment->Fill(9);
01422 
01423   return true;
01424 }
01425 //
01426 void CSCEfficiency::returnTypes(CSCDetId & id, int &ec, int &st, int &rg, int &ch, int &secondRing){
01427   ec = id.endcap()-1;
01428   st = id.station()-1;
01429   rg = id.ring()-1;
01430   secondRing = -1;
01431   if(1==id.station() && (4==id.ring() || 1==id.ring()) ){
01432     rg = 0;
01433     secondRing = 3;
01434   }
01435   ch = id.chamber()-FirstCh;
01436 }
01437 
01438 //
01439 void CSCEfficiency::getFromFTS( const FreeTrajectoryState& fts,
01440                                 CLHEP::Hep3Vector& p3, CLHEP::Hep3Vector& r3,
01441                                 int& charge, AlgebraicSymMatrix66& cov){
01442   
01443   GlobalVector p3GV = fts.momentum();
01444   GlobalPoint r3GP = fts.position();
01445 
01446   p3.set(p3GV.x(), p3GV.y(), p3GV.z());
01447   r3.set(r3GP.x(), r3GP.y(), r3GP.z());
01448 
01449   charge = fts.charge();
01450   cov = fts.hasError() ? fts.cartesianError().matrix() : AlgebraicSymMatrix66();
01451 
01452 }
01453 
01454 FreeTrajectoryState CSCEfficiency::getFromCLHEP(const CLHEP::Hep3Vector& p3, const CLHEP::Hep3Vector& r3,
01455                                                        int charge, const AlgebraicSymMatrix66& cov,
01456                                                        const MagneticField* field){
01457 
01458   GlobalVector p3GV(p3.x(), p3.y(), p3.z());
01459   GlobalPoint r3GP(r3.x(), r3.y(), r3.z());
01460   GlobalTrajectoryParameters tPars(r3GP, p3GV, charge, field);
01461 
01462   CartesianTrajectoryError tCov(cov);
01463 
01464   return cov.kRows == 6 ? FreeTrajectoryState(tPars, tCov) : FreeTrajectoryState(tPars) ;
01465 }
01466 
01467 void CSCEfficiency::linearExtrapolation(GlobalPoint initialPosition ,GlobalVector initialDirection, 
01468                                         float zSurface, std::vector <float> &posZY){
01469   double paramLine =  lineParameter(initialPosition.z(), zSurface, initialDirection.z());
01470   double xPosition = extrapolate1D(initialPosition.x(), initialDirection.x(),paramLine);
01471   double yPosition = extrapolate1D(initialPosition.y(), initialDirection.y(),paramLine);
01472   posZY.clear();
01473   posZY.push_back(xPosition);
01474   posZY.push_back(yPosition);
01475 }
01476 //
01477 double CSCEfficiency::extrapolate1D(double initPosition, double initDirection, double parameterOfTheLine){
01478   double extrapolatedPosition = initPosition + initDirection*parameterOfTheLine;
01479   return extrapolatedPosition;
01480 }
01481 //
01482 double CSCEfficiency::lineParameter(double initZPosition, double destZPosition, double initZDirection){
01483   double paramLine = (destZPosition-initZPosition)/initZDirection;
01484   return paramLine;
01485 }
01486 //
01487 void CSCEfficiency::chooseDirection(CLHEP::Hep3Vector & innerPosition, CLHEP::Hep3Vector & outerPosition){
01488 
01489   //---- Be careful with trigger conditions too
01490   if(!isIPdata){
01491     float dy = outerPosition.y() - innerPosition.y();
01492     float dz = outerPosition.z() - innerPosition.z();
01493     if(isBeamdata){
01494       if(dz>0){
01495         alongZ = true;
01496       }
01497       else{
01498         alongZ = false;
01499       }
01500     }
01501     else{//cosmics
01502       if(dy/dz>0){
01503         alongZ = false;
01504       }
01505       else{
01506         alongZ = true;
01507       }
01508     }
01509   }
01510 }
01511 //
01512 const Propagator* CSCEfficiency::propagator(std::string propagatorName) const {
01513   return &*theService->propagator(propagatorName);
01514 }
01515   
01516 //
01517 TrajectoryStateOnSurface CSCEfficiency::propagate(FreeTrajectoryState & ftsStart, const BoundPlane &bpDest){
01518   TrajectoryStateOnSurface tSOSDest;
01519   std::string propagatorName;
01520 /*
01521 // it would work if cosmic muons had properly assigned direction...
01522   bool dzPositive = bpDest.position().z() - ftsStart.position().z() > 0 ? true : false;
01523  //---- Be careful with trigger conditions too
01524   if(!isIPdata){
01525     bool rightDirection = !(alongZ^dzPositive);
01526     if(rightDirection){
01527       if(printalot) std::cout<<" propagate along momentum"<<std::endl;
01528       propagatorName = "SteppingHelixPropagatorAlong";
01529     }
01530     else{
01531       if(printalot) std::cout<<" propagate opposite momentum"<<std::endl;
01532       propagatorName = "SteppingHelixPropagatorOpposite";
01533     }
01534   }
01535   else{
01536     if(printalot) std::cout<<" propagate any (momentum)"<<std::endl;
01537     propagatorName = "SteppingHelixPropagatorAny";
01538   }
01539 */
01540   propagatorName = "SteppingHelixPropagatorAny";
01541   tSOSDest = propagator(propagatorName)->propagate(ftsStart, bpDest);   
01542   return tSOSDest;
01543 }
01544 //
01545 bool CSCEfficiency::applyTrigger(edm::Handle<edm::TriggerResults> &hltR,
01546                                  const edm::TriggerNames & triggerNames){
01547   bool triggerPassed = true;
01548   std::vector<std::string>  hlNames=triggerNames.triggerNames();
01549   pointToTriggers.clear();
01550   for(size_t imyT = 0;imyT<myTriggers.size();++imyT){
01551     for (size_t iT=0; iT<hlNames.size(); ++iT) {
01552       //std::cout<<" iT = "<<iT<<" hlNames[iT] = "<<hlNames[iT]<<
01553       //" : wasrun = "<<hltR->wasrun(iT)<<" accept = "<<
01554       //         hltR->accept(iT)<<" !error = "<< 
01555       //        !hltR->error(iT)<<std::endl; 
01556       if(!imyT){       
01557         if(hltR->wasrun(iT) &&
01558            hltR->accept(iT) &&
01559            !hltR->error(iT) ){
01560            TriggersFired->Fill(iT); 
01561         }
01562       }
01563       if(hlNames[iT]==myTriggers[imyT]){
01564         pointToTriggers.push_back(iT);
01565         if(imyT){
01566           break;
01567         }
01568       }
01569     }
01570   }
01571   if(pointToTriggers.size()!=myTriggers.size()){
01572     pointToTriggers.clear();
01573     if(printalot){
01574       std::cout<<" Not all trigger names found - all trigger specifications will be ignored. Check your cfg file!"<<std::endl;
01575     }
01576   }
01577   else{
01578     if(pointToTriggers.size()){
01579       if(printalot){
01580         std::cout<<"The following triggers will be required in the event: "<<std::endl;
01581         for(size_t imyT =0; imyT <pointToTriggers.size();++imyT){
01582           std::cout<<"  "<<hlNames[pointToTriggers[imyT]];
01583         }
01584         std::cout<<std::endl;
01585         std::cout<<" in condition (AND/OR) : "<<!andOr<<"/"<<andOr<<std::endl;
01586       }
01587     }
01588   }
01589 
01590   if (hltR.isValid()) {
01591     if(!pointToTriggers.size()){
01592       if(printalot){
01593         std::cout<<" No triggers specified in the configuration or all ignored - no trigger information will be considered"<<std::endl;
01594       }
01595     }
01596     for(size_t imyT =0; imyT <pointToTriggers.size();++imyT){
01597       if(hltR->wasrun(pointToTriggers[imyT]) && 
01598          hltR->accept(pointToTriggers[imyT]) && 
01599          !hltR->error(pointToTriggers[imyT]) ){
01600         triggerPassed = true;
01601         if(andOr){
01602           break;
01603         }
01604       }
01605       else{
01606         triggerPassed = false;
01607         if(!andOr){
01608           triggerPassed = false;
01609           break;
01610         }
01611       }
01612     }
01613   }
01614   else{
01615     if(printalot){
01616       std::cout<<" TriggerResults handle returns invalid state?! No trigger information will be considered"<<std::endl;
01617     }
01618   }
01619   if(printalot){
01620     std::cout<<" Trigger passed: "<<triggerPassed<<std::endl;
01621   }
01622   return triggerPassed;
01623 }
01624 //
01625 
01626 // Constructor
01627 CSCEfficiency::CSCEfficiency(const ParameterSet& pset){
01628 
01629   // const float Xmin = -70;
01630   //const float Xmax = 70;
01631   //const int nXbins = int(4.*(Xmax - Xmin));
01632   const float Ymin = -165;
01633   const float Ymax = 165;
01634   const int nYbins = int((Ymax - Ymin)/2);
01635   const float Layer_min = -0.5;
01636   const float Layer_max = 9.5;
01637   const int nLayer_bins = int(Layer_max - Layer_min);
01638   //
01639 
01640   //---- Get the input parameters
01641   printout_NEvents  = pset.getUntrackedParameter<unsigned int>("printout_NEvents",0);
01642   rootFileName     = pset.getUntrackedParameter<string>("rootFileName","cscHists.root");
01643 
01644   isData  = pset.getUntrackedParameter<bool>("runOnData",true);// 
01645   isIPdata  = pset.getUntrackedParameter<bool>("IPdata",false);// 
01646   isBeamdata  = pset.getUntrackedParameter<bool>("Beamdata",false);//
01647   getAbsoluteEfficiency  = pset.getUntrackedParameter<bool>("getAbsoluteEfficiency",true);//
01648   useDigis = pset.getUntrackedParameter<bool>("useDigis", true);// 
01649   distanceFromDeadZone = pset.getUntrackedParameter<double>("distanceFromDeadZone", 10.);// 
01650   minP = pset.getUntrackedParameter<double>("minP",20.);//
01651   maxP = pset.getUntrackedParameter<double>("maxP",100.);//
01652   maxNormChi2 = pset.getUntrackedParameter<double>("maxNormChi2", 3.);//
01653   minTrackHits = pset.getUntrackedParameter<unsigned int>("minTrackHits",10);//
01654 
01655   applyIPangleCuts = pset.getUntrackedParameter<bool>("applyIPangleCuts", false);// 
01656     local_DY_DZ_Max = pset.getUntrackedParameter<double>("local_DY_DZ_Max",-0.1);//
01657       local_DY_DZ_Min = pset.getUntrackedParameter<double>("local_DY_DZ_Min",-0.8);//
01658         local_DX_DZ_Max = pset.getUntrackedParameter<double>("local_DX_DZ_Max",0.2);//
01659 
01660   alctDigiTag_  = pset.getParameter<edm::InputTag>("alctDigiTag") ;
01661   clctDigiTag_  = pset.getParameter<edm::InputTag>("clctDigiTag") ;
01662   corrlctDigiTag_  = pset.getParameter<edm::InputTag>("corrlctDigiTag") ;
01663   stripDigiTag_  = pset.getParameter<edm::InputTag>("stripDigiTag") ;
01664   wireDigiTag_ = pset.getParameter<edm::InputTag>("wireDigiTag") ;
01665   rechitDigiTag_ = pset.getParameter<edm::InputTag>("rechitDigiTag") ;
01666   segmentDigiTag_ = pset.getParameter<edm::InputTag>("segmentDigiTag") ;
01667   simHitTag     = pset.getParameter<edm::InputTag>("simHitTag");
01668   tracksTag = pset.getParameter< edm::InputTag >("tracksTag");
01669 
01670   ParameterSet serviceParameters = pset.getParameter<ParameterSet>("ServiceParameters");
01671   // maybe use the service for getting magnetic field, propagators, etc. ...
01672   theService        = new MuonServiceProxy(serviceParameters);
01673 
01674   // Trigger
01675   useTrigger =  pset.getUntrackedParameter<bool>("useTrigger", false);
01676   hlTriggerResults_ = pset.getParameter<edm::InputTag> ("HLTriggerResults");
01677   myTriggers = pset.getParameter<std::vector <std::string> >("myTriggers");
01678   andOr =  pset.getUntrackedParameter<bool>("andOr");
01679   pointToTriggers.clear();
01680 
01681 
01682   //---- set counter to zero
01683   nEventsAnalyzed = 0;
01684   //---- set presence of magnetic field
01685   magField = true;
01686   //
01687   std::string Path = "AllChambers/";
01688   std::string FullName;
01689   //---- File with output histograms 
01690   theFile = new TFile(rootFileName.c_str(), "RECREATE");
01691   theFile->cd();
01692   //---- Book histograms for the analysis
01693   char SpecName[50];
01694   
01695   sprintf(SpecName,"DataFlow"); 
01696   DataFlow =  
01697     new TH1F(SpecName,"Data flow;condition number;entries",40,-0.5,39.5);
01698   //
01699   sprintf(SpecName,"TriggersFired"); 
01700   TriggersFired =
01701     new TH1F(SpecName,"Triggers fired;trigger number;entries",140,-0.5,139.5);
01702   //
01703   int Chan = 50;
01704   float minChan = -0.5;
01705   float maxChan = 49.5;
01706   //
01707   sprintf(SpecName,"ALCTPerEvent");    
01708   ALCTPerEvent = new TH1F(SpecName,"ALCTs per event;N digis;entries",Chan,minChan,maxChan);
01709   //
01710   sprintf(SpecName,"CLCTPerEvent");    
01711   CLCTPerEvent = new TH1F(SpecName,"CLCTs per event;N digis;entries",Chan,minChan,maxChan);
01712   //
01713     sprintf(SpecName,"recHitsPerEvent");    
01714   recHitsPerEvent = new TH1F(SpecName,"RecHits per event;N digis;entries",150,-0.5,149.5);
01715   //
01716   sprintf(SpecName,"segmentsPerEvent");    
01717   segmentsPerEvent = new TH1F(SpecName,"segments per event;N digis;entries",Chan,minChan,maxChan);
01718   //
01719   //---- Book groups of histograms (for any chamber)
01720 
01721   map<std::string,bool>::iterator iter;  
01722   for(int ec = 0;ec<2;++ec){
01723     for(int st = 0;st<4;++st){
01724       theFile->cd();
01725       sprintf(SpecName,"Stations__E%d_S%d",ec+1, st+1);
01726       theFile->mkdir(SpecName);
01727       theFile->cd(SpecName);
01728 
01729       //
01730       sprintf(SpecName,"segmentChi2_ndf_St%d",st+1);
01731       StHist[ec][st].segmentChi2_ndf = 
01732         new TH1F(SpecName,"Chi2/ndf of a segment;chi2/ndf;entries",100,0.,20.);
01733       //
01734       sprintf(SpecName,"hitsInSegment_St%d",st+1);
01735       StHist[ec][st].hitsInSegment = 
01736         new TH1F(SpecName,"Number of hits in a segment;nHits;entries",7,-0.5,6.5);
01737       //
01738       Chan = 170;
01739       minChan = 0.85;
01740       maxChan = 2.55;
01741       //
01742       sprintf(SpecName,"AllSegments_eta_St%d",st+1);
01743       StHist[ec][st].AllSegments_eta = 
01744         new TH1F(SpecName,"All segments in eta;eta;entries",Chan,minChan,maxChan);
01745       //
01746       sprintf(SpecName,"EfficientSegments_eta_St%d",st+1);
01747       StHist[ec][st].EfficientSegments_eta = 
01748         new TH1F(SpecName,"Efficient segments in eta;eta;entries",Chan,minChan,maxChan);
01749       //
01750       sprintf(SpecName,"ResidualSegments_St%d",st+1);
01751       StHist[ec][st].ResidualSegments = 
01752         new TH1F(SpecName,"Residual (segments);residual,cm;entries",75,0.,15.);
01753       //
01754       Chan = 200;
01755       minChan = -800.;
01756       maxChan = 800.;
01757       int Chan2 = 200;
01758       float minChan2 = -800.;
01759       float maxChan2 = 800.;
01760       
01761       sprintf(SpecName,"EfficientSegments_XY_St%d",st+1);
01762       StHist[ec][st].EfficientSegments_XY = new TH2F(SpecName,"Efficient segments in XY;X;Y",
01763                                                      Chan,minChan,maxChan,Chan2,minChan2,maxChan2);
01764       sprintf(SpecName,"InefficientSegments_XY_St%d",st+1);
01765       StHist[ec][st].InefficientSegments_XY = new TH2F(SpecName,"Inefficient segments in XY;X;Y",
01766                                                        Chan,minChan,maxChan,Chan2,minChan2,maxChan2);
01767       //
01768       Chan = 80;
01769       minChan = 0;
01770       maxChan = 3.2;
01771       sprintf(SpecName,"EfficientALCT_momTheta_St%d",st+1);
01772       StHist[ec][st].EfficientALCT_momTheta = new TH1F(SpecName,"Efficient ALCT in theta (momentum);theta, rad;entries",
01773                                                        Chan,minChan,maxChan);
01774       //
01775       sprintf(SpecName,"InefficientALCT_momTheta_St%d",st+1);
01776       StHist[ec][st].InefficientALCT_momTheta = new TH1F(SpecName,"Inefficient ALCT in theta (momentum);theta, rad;entries",
01777                                                          Chan,minChan,maxChan);
01778       //
01779       Chan = 160;
01780       minChan = -3.2;
01781       maxChan = 3.2;
01782       sprintf(SpecName,"EfficientCLCT_momPhi_St%d",st+1);
01783       StHist[ec][st].EfficientCLCT_momPhi = new TH1F(SpecName,"Efficient CLCT in phi (momentum);phi, rad;entries",
01784                                                      Chan,minChan,maxChan);
01785       //
01786       sprintf(SpecName,"InefficientCLCT_momPhi_St%d",st+1);
01787       StHist[ec][st].InefficientCLCT_momPhi = new TH1F(SpecName,"Inefficient CLCT in phi (momentum);phi, rad;entries",
01788                                                        Chan,minChan,maxChan);
01789       //
01790       theFile->cd();
01791       for(int rg = 0;rg<3;++rg){
01792         if(0!=st && rg>1){
01793           continue;
01794         }
01795         else if(1==rg && 3==st){
01796           continue;
01797         }
01798         for(int iChamber=FirstCh;iChamber<FirstCh+NumCh;iChamber++){
01799           if(0!=st && 0==rg && iChamber >18){
01800             continue;
01801           }
01802           theFile->cd();
01803           sprintf(SpecName,"Chambers__E%d_S%d_R%d_Chamber_%d",ec+1, st+1, rg+1,iChamber);
01804           theFile->mkdir(SpecName);
01805           theFile->cd(SpecName);
01806           //
01807 
01808           sprintf(SpecName,"EfficientRechits_inSegment_Ch%d",iChamber);
01809           ChHist[ec][st][rg][iChamber-FirstCh].EfficientRechits_inSegment = 
01810             new TH1F(SpecName,"Existing RecHit given a segment;layers (1-6);entries",nLayer_bins,Layer_min,Layer_max);
01811           //
01812           sprintf(SpecName,"InefficientSingleHits_Ch%d",iChamber);
01813           ChHist[ec][st][rg][iChamber-FirstCh].InefficientSingleHits = 
01814             new TH1F(SpecName,"Single RecHits not in the segment;layers (1-6);entries ",nLayer_bins,Layer_min,Layer_max);
01815           //
01816           sprintf(SpecName,"AllSingleHits_Ch%d",iChamber);
01817           ChHist[ec][st][rg][iChamber-FirstCh].AllSingleHits = 
01818             new TH1F(SpecName,"Single RecHits given a segment; layers (1-6);entries",nLayer_bins,Layer_min,Layer_max);
01819           //
01820           sprintf(SpecName,"digiAppearanceCount_Ch%d",iChamber);
01821           ChHist[ec][st][rg][iChamber-FirstCh].digiAppearanceCount = 
01822             new TH1F(SpecName,"Digi appearance (no-yes): segment(0,1), ALCT(2,3), CLCT(4,5), CorrLCT(6,7); digi type;entries",
01823                      8,-0.5,7.5);
01824           //
01825           Chan = 100;
01826           minChan = -1.1;
01827           maxChan = 0.9;
01828           sprintf(SpecName,"EfficientALCT_dydz_Ch%d",iChamber);
01829           ChHist[ec][st][rg][iChamber-FirstCh].EfficientALCT_dydz = 
01830             new TH1F(SpecName,"Efficient ALCT; local dy/dz (ME 3 and 4 flipped);entries",
01831                      Chan, minChan, maxChan);
01832           //
01833           sprintf(SpecName,"InefficientALCT_dydz_Ch%d",iChamber);
01834           ChHist[ec][st][rg][iChamber-FirstCh].InefficientALCT_dydz = 
01835             new TH1F(SpecName,"Inefficient ALCT; local dy/dz (ME 3 and 4 flipped);entries",
01836                      Chan, minChan, maxChan);
01837           //
01838           Chan = 100;
01839           minChan = -1.;
01840           maxChan = 1.0;
01841           sprintf(SpecName,"EfficientCLCT_dxdz_Ch%d",iChamber);
01842           ChHist[ec][st][rg][iChamber-FirstCh].EfficientCLCT_dxdz = 
01843             new TH1F(SpecName,"Efficient CLCT; local dxdz;entries",
01844                      Chan, minChan, maxChan);
01845           //
01846           sprintf(SpecName,"InefficientCLCT_dxdz_Ch%d",iChamber);
01847           ChHist[ec][st][rg][iChamber-FirstCh].InefficientCLCT_dxdz = 
01848             new TH1F(SpecName,"Inefficient CLCT; local dxdz;entries",
01849                      Chan, minChan, maxChan);
01850           //
01851           sprintf(SpecName,"EfficientRechits_good_Ch%d",iChamber);
01852           ChHist[ec][st][rg][iChamber-FirstCh].EfficientRechits_good = 
01853             new TH1F(SpecName,"Existing RecHit - sensitive area only;layers (1-6);entries",nLayer_bins,Layer_min,Layer_max);
01854           //
01855           sprintf(SpecName,"EfficientStrips_Ch%d",iChamber);
01856           ChHist[ec][st][rg][iChamber-FirstCh].EfficientStrips = 
01857             new TH1F(SpecName,"Existing strip;layer (1-6); entries",nLayer_bins,Layer_min,Layer_max);
01858           //
01859           sprintf(SpecName,"EfficientWireGroups_Ch%d",iChamber);
01860           ChHist[ec][st][rg][iChamber-FirstCh].EfficientWireGroups = 
01861             new TH1F(SpecName,"Existing WireGroups;layer (1-6); entries ",nLayer_bins,Layer_min,Layer_max);
01862           //
01863           sprintf(SpecName,"StripWiresCorrelations_Ch%d",iChamber);
01864           ChHist[ec][st][rg][iChamber-FirstCh].StripWiresCorrelations = 
01865             new TH1F(SpecName,"StripWire correlations;; entries ",5,0.5,5.5);
01866           //
01867           Chan = 80;
01868           minChan = 0;
01869           maxChan = 3.2;
01870           sprintf(SpecName,"NoWires_momTheta_Ch%d",iChamber);
01871           ChHist[ec][st][rg][iChamber-FirstCh].NoWires_momTheta = 
01872             new TH1F(SpecName,"No wires (all strips present) - in theta (momentum);theta, rad;entries",
01873                      Chan,minChan,maxChan);
01874           //
01875           Chan = 160;
01876           minChan = -3.2;
01877           maxChan = 3.2;
01878           sprintf(SpecName,"NoStrips_momPhi_Ch%d",iChamber);
01879           ChHist[ec][st][rg][iChamber-FirstCh].NoStrips_momPhi = 
01880             new TH1F(SpecName,"No strips (all wires present) - in phi (momentum);phi, rad;entries",
01881                      Chan,minChan,maxChan);
01882           //
01883           for(int iLayer=0; iLayer<6;iLayer++){
01884             sprintf(SpecName,"Y_InefficientRecHits_inSegment_Ch%d_L%d",iChamber,iLayer);
01885             ChHist[ec][st][rg][iChamber-FirstCh].Y_InefficientRecHits_inSegment.push_back
01886               (new TH1F(SpecName,"Missing RecHit/layer in a segment (local system, whole chamber);Y, cm; entries",
01887                         nYbins,Ymin, Ymax));
01888             //
01889             sprintf(SpecName,"Y_EfficientRecHits_inSegment_Ch%d_L%d",iChamber,iLayer);
01890             ChHist[ec][st][rg][iChamber-FirstCh].Y_EfficientRecHits_inSegment.push_back
01891               (new TH1F(SpecName,"Efficient (extrapolated from the segment) RecHit/layer in a segment (local system, whole chamber);Y, cm; entries",
01892                         nYbins,Ymin, Ymax));
01893             //
01894             Chan = 200;
01895             minChan = -0.2;
01896             maxChan = 0.2;
01897             sprintf(SpecName,"Phi_InefficientRecHits_inSegment_Ch%d_L%d",iChamber,iLayer);
01898             ChHist[ec][st][rg][iChamber-FirstCh].Phi_InefficientRecHits_inSegment.push_back
01899               (new TH1F(SpecName,"Missing RecHit/layer in a segment (local system, whole chamber);Phi, rad; entries",
01900                         Chan, minChan, maxChan));
01901             //
01902             sprintf(SpecName,"Phi_EfficientRecHits_inSegment_Ch%d_L%d",iChamber,iLayer);
01903             ChHist[ec][st][rg][iChamber-FirstCh].Phi_EfficientRecHits_inSegment.push_back
01904               (new TH1F(SpecName,"Efficient (extrapolated from the segment) in a segment (local system, whole chamber);Phi, rad; entries",
01905                         Chan, minChan, maxChan));
01906             
01907           }
01908           //
01909           sprintf(SpecName,"Sim_Rechits_Ch%d",iChamber);
01910           ChHist[ec][st][rg][iChamber-FirstCh].SimRechits = 
01911             new TH1F(SpecName,"Existing RecHit (Sim);layers (1-6);entries",nLayer_bins,Layer_min,Layer_max);
01912           //
01913           sprintf(SpecName,"Sim_Simhits_Ch%d",iChamber);
01914           ChHist[ec][st][rg][iChamber-FirstCh].SimSimhits = 
01915             new TH1F(SpecName,"Existing SimHit (Sim);layers (1-6);entries",nLayer_bins,Layer_min,Layer_max);
01916           //
01917           /*
01918           sprintf(SpecName,"Sim_Rechits_each_Ch%d",iChamber);
01919           ChHist[ec][st][rg][iChamber-FirstCh].SimRechits_each = 
01920             new TH1F(SpecName,"Existing RecHit (Sim), each;layers (1-6);entries",nLayer_bins,Layer_min,Layer_max);
01921           //
01922           sprintf(SpecName,"Sim_Simhits_each_Ch%d",iChamber);
01923           ChHist[ec][st][rg][iChamber-FirstCh].SimSimhits_each = 
01924             new TH1F(SpecName,"Existing SimHit (Sim), each;layers (1-6);entries",nLayer_bins,Layer_min,Layer_max);
01925           */
01926           theFile->cd();
01927         }
01928       }
01929     }
01930   }
01931 }
01932 
01933 // Destructor
01934 CSCEfficiency::~CSCEfficiency(){
01935   if (theService) delete theService; 
01936   // Write the histos to a file
01937   theFile->cd();
01938   //
01939   char SpecName[20];
01940   std::vector<float> bins, Efficiency, EffError;
01941   std::vector<float> eff(2);
01942 
01943   //---- loop over chambers
01944   std::map <std::string, bool> chamberTypes;
01945   chamberTypes["ME11"] = false;
01946   chamberTypes["ME12"] = false;
01947   chamberTypes["ME13"] = false;
01948   chamberTypes["ME21"] = false;
01949   chamberTypes["ME22"] = false;
01950   chamberTypes["ME31"] = false;
01951   chamberTypes["ME32"] = false;
01952   chamberTypes["ME41"] = false;
01953   
01954   map<std::string,bool>::iterator iter;  
01955   std::cout<<" Writing proper histogram structure (patience)..."<<std::endl;
01956   for(int ec = 0;ec<2;++ec){
01957     for(int st = 0;st<4;++st){
01958       sprintf(SpecName,"Stations__E%d_S%d",ec+1, st+1);
01959       theFile->cd(SpecName);
01960       StHist[ec][st].segmentChi2_ndf->Write();
01961       StHist[ec][st].hitsInSegment->Write();
01962       StHist[ec][st].AllSegments_eta->Write();
01963       StHist[ec][st].EfficientSegments_eta->Write();
01964       StHist[ec][st].ResidualSegments->Write();
01965       StHist[ec][st].EfficientSegments_XY->Write();
01966       StHist[ec][st].InefficientSegments_XY->Write();
01967       StHist[ec][st].EfficientALCT_momTheta->Write();
01968       StHist[ec][st].InefficientALCT_momTheta->Write();
01969       StHist[ec][st].EfficientCLCT_momPhi->Write();
01970       StHist[ec][st].InefficientCLCT_momPhi->Write();
01971       for(int rg = 0;rg<3;++rg){
01972         if(0!=st && rg>1){
01973           continue;
01974         }
01975         else if(1==rg && 3==st){
01976           continue;
01977         }
01978         for(int iChamber=FirstCh;iChamber<FirstCh+NumCh;iChamber++){
01979           if(0!=st && 0==rg && iChamber >18){
01980             continue;
01981           }
01982           sprintf(SpecName,"Chambers__E%d_S%d_R%d_Chamber_%d",ec+1, st+1, rg+1,iChamber);
01983           theFile->cd(SpecName);
01984           
01985           ChHist[ec][st][rg][iChamber-FirstCh].EfficientRechits_inSegment->Write();
01986           ChHist[ec][st][rg][iChamber-FirstCh].AllSingleHits->Write();
01987           ChHist[ec][st][rg][iChamber-FirstCh].digiAppearanceCount->Write();
01988           ChHist[ec][st][rg][iChamber-FirstCh].EfficientALCT_dydz->Write();
01989           ChHist[ec][st][rg][iChamber-FirstCh].InefficientALCT_dydz->Write();
01990           ChHist[ec][st][rg][iChamber-FirstCh].EfficientCLCT_dxdz->Write();
01991           ChHist[ec][st][rg][iChamber-FirstCh].InefficientCLCT_dxdz->Write();
01992           ChHist[ec][st][rg][iChamber-FirstCh].InefficientSingleHits->Write();
01993           ChHist[ec][st][rg][iChamber-FirstCh].EfficientRechits_good->Write();
01994           ChHist[ec][st][rg][iChamber-FirstCh].EfficientStrips->Write();
01995           ChHist[ec][st][rg][iChamber-FirstCh].StripWiresCorrelations->Write();
01996           ChHist[ec][st][rg][iChamber-FirstCh].NoWires_momTheta->Write();
01997           ChHist[ec][st][rg][iChamber-FirstCh].NoStrips_momPhi->Write();
01998           ChHist[ec][st][rg][iChamber-FirstCh].EfficientWireGroups->Write();
01999           for(unsigned int iLayer = 0; iLayer< 6; iLayer++){
02000             ChHist[ec][st][rg][iChamber-FirstCh].Y_InefficientRecHits_inSegment[iLayer]->Write();
02001             ChHist[ec][st][rg][iChamber-FirstCh].Y_EfficientRecHits_inSegment[iLayer]->Write();
02002             ChHist[ec][st][rg][iChamber-FirstCh].Phi_InefficientRecHits_inSegment[iLayer]->Write();
02003             ChHist[ec][st][rg][iChamber-FirstCh].Phi_EfficientRecHits_inSegment[iLayer]->Write();
02004           }
02005           ChHist[ec][st][rg][iChamber-FirstCh].SimRechits->Write();
02006           ChHist[ec][st][rg][iChamber-FirstCh].SimSimhits->Write();
02007           /*
02008           ChHist[ec][st][rg][iChamber-FirstCh].SimRechits_each->Write();
02009           ChHist[ec][st][rg][iChamber-FirstCh].SimSimhits_each->Write();
02010           */
02011           //
02012           theFile->cd(SpecName);
02013           theFile->cd();
02014         }
02015       }
02016     }
02017   }
02018   //
02019   sprintf(SpecName,"AllChambers");
02020   theFile->mkdir(SpecName);
02021   theFile->cd(SpecName);
02022   DataFlow->Write(); 
02023   TriggersFired->Write();
02024   ALCTPerEvent->Write();
02025   CLCTPerEvent->Write();
02026   recHitsPerEvent->Write();
02027   segmentsPerEvent->Write();
02028   //
02029   theFile->cd(SpecName);
02030   //---- Close the file
02031   theFile->Close();
02032 }
02033 
02034 // ------------ method called once each job just before starting event loop  ------------
02035 void
02036 CSCEfficiency::beginJob()
02037 {
02038 }
02039 
02040 // ------------ method called once each job just after ending the event loop  ------------
02041 void
02042 CSCEfficiency::endJob() {
02043 }
02044 
02045 DEFINE_FWK_MODULE(CSCEfficiency);
02046
Generated for CMSSW Reference Manual by doxygen 1.7.3