CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_0/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 
00788     const std::vector<float> LayerBounds = layerGeom->parameters ();
00789     std::vector<CSCWireDigi>::const_iterator digiItr = (*j).second.first;
00790     std::vector<CSCWireDigi>::const_iterator last = (*j).second.second;
00791     //
00792     for( ; digiItr != last; ++digiItr) {
00793       std::pair < int, float > WG_pos(digiItr->getWireGroup(), layerGeom->yOfWireGroup(digiItr->getWireGroup())); 
00794       std::pair <std::pair < int, float >, int >  LayerSignal(WG_pos, digiItr->getTimeBin()); 
00795       
00796       //---- AllWG contains basic information about WG (WG number and Y-position, time bin)
00797       allWG[id.endcap()-1][id.station()-1][id.ring()-1][id.chamber()-FirstCh]
00798         [id.layer()-1].push_back(LayerSignal);
00799       if(printalot){
00800         //std::cout<<" WG check : "<<std::endl;
00801         //printf("\t\tendcap/station/ring/chamber/layer: %i/%i/%i/%i/%i\n",id.endcap(),id.station(),id.ring(),id.chamber(),id.layer());
00802         //std::cout<<" WG size = "<<allWG[id.endcap()-1][id.station()-1][id.ring()-1][id.chamber()-FirstCh]
00803         //[id.layer()-1].size()<<std::endl;
00804       }
00805     }
00806   }
00807 }
00808 void CSCEfficiency::fillStrips_info(edm::Handle<CSCStripDigiCollection> &strips){
00809   //---- STRIPS
00810   for (CSCStripDigiCollection::DigiRangeIterator j=strips->begin(); j!=strips->end(); j++) {
00811     CSCDetId id = (CSCDetId)(*j).first;
00812     int largestADCValue = -1;
00813     std::vector<CSCStripDigi>::const_iterator digiItr = (*j).second.first;
00814     std::vector<CSCStripDigi>::const_iterator last = (*j).second.second;
00815     for( ; digiItr != last; ++digiItr) {
00816       int maxADC=largestADCValue;
00817       int myStrip = digiItr->getStrip();
00818       std::vector<int> myADCVals = digiItr->getADCCounts();
00819       float thisPedestal = 0.5*(float)(myADCVals[0]+myADCVals[1]);
00820       float threshold = 13.3 ;
00821       float diff = 0.;
00822       float peakADC  = -1000.;
00823       for (unsigned int iCount = 0; iCount < myADCVals.size(); iCount++) {
00824         diff = (float)myADCVals[iCount]-thisPedestal;
00825         if (diff > threshold) { 
00826           if (myADCVals[iCount] > largestADCValue) {
00827             largestADCValue = myADCVals[iCount];
00828           }
00829         }
00830         if (diff > threshold && diff > peakADC) {
00831           peakADC  = diff;
00832         }
00833       }
00834       if(largestADCValue>maxADC){// FIX IT!!!
00835         maxADC = largestADCValue;
00836         std::pair <int, float> LayerSignal (myStrip, peakADC);
00837         
00838         //---- AllStrips contains basic information about strips 
00839         //---- (strip number and peak signal for most significant strip in the layer) 
00840         allStrips[id.endcap()-1][id.station()-1][id.ring()-1][id.chamber()-1][id.layer()-1].clear();
00841         allStrips[id.endcap()-1][id.station()-1][id.ring()-1][id.chamber()-1][id.layer()-1].push_back(LayerSignal);
00842       }
00843     }
00844   }
00845 }
00846 void CSCEfficiency::fillSimhit_info(edm::Handle<edm::PSimHitContainer> &simhits){
00847   //---- SIMHITS
00848   edm::PSimHitContainer::const_iterator dSHsimIter;
00849   for (dSHsimIter = simhits->begin(); dSHsimIter != simhits->end(); dSHsimIter++){
00850     // Get DetID for this simHit:
00851     CSCDetId sId = (CSCDetId)(*dSHsimIter).detUnitId();
00852     std::pair <LocalPoint, int> simHitPos((*dSHsimIter).localPosition(), (*dSHsimIter).particleType());  
00853     allSimhits[sId.endcap()-1][sId.station()-1][sId.ring()-1][sId.chamber()-FirstCh][sId.layer()-1].push_back(simHitPos);
00854   }
00855 }
00856 //
00857 void CSCEfficiency::fillRechitsSegments_info(edm::Handle<CSCRecHit2DCollection> &rechits,
00858                                              edm::Handle<CSCSegmentCollection> &segments,
00859                                              edm::ESHandle<CSCGeometry> &cscGeom
00860                                              ){
00861   //---- RECHITS AND SEGMENTS
00862   //---- Loop over rechits 
00863   if (printalot){ 
00864     //printf("\tGet the recHits collection.\t ");
00865     printf("  The size of the rechit collection is %i\n",int(rechits->size()));
00866     //printf("\t...start loop over rechits...\n");
00867   }
00868   recHitsPerEvent->Fill(rechits->size());
00869   //---- Build iterator for rechits and loop :
00870   CSCRecHit2DCollection::const_iterator recIt;
00871   for (recIt = rechits->begin(); recIt != rechits->end(); recIt++) {
00872     //---- Find chamber with rechits in CSC 
00873     CSCDetId id = (CSCDetId)(*recIt).cscDetId();
00874     if (printalot){
00875       const CSCLayer* csclayer = cscGeom->layer( id);
00876       LocalPoint rhitlocal = (*recIt).localPosition();  
00877       LocalError rerrlocal = (*recIt).localPositionError();  
00878       GlobalPoint rhitglobal= csclayer->toGlobal(rhitlocal);
00879       printf("\t\tendcap/station/ring/chamber/layer: %i/%i/%i/%i/%i\n",id.endcap(),id.station(),id.ring(),id.chamber(),id.layer());
00880       printf("\t\tx,y,z: %f, %f, %f\texx,eey,exy: %f, %f, %f\tglobal x,y,z: %f, %f, %f \n",
00881              rhitlocal.x(), rhitlocal.y(), rhitlocal.z(), rerrlocal.xx(), rerrlocal.yy(), rerrlocal.xy(),
00882              rhitglobal.x(), rhitglobal.y(), rhitglobal.z());
00883     }
00884     std::pair <LocalPoint, bool> recHitPos((*recIt).localPosition(), false);  
00885     allRechits[id.endcap()-1][id.station()-1][id.ring()-1][id.chamber()-FirstCh][id.layer()-1].push_back(recHitPos);
00886   }
00887   //---- "Empty" chambers
00888   for(int iE=0;iE<2;iE++){
00889     for(int iS=0;iS<4;iS++){
00890       for(int iR=0;iR<4;iR++){
00891         for(int iC=0;iC<NumCh;iC++){
00892           int numLayers = 0;
00893           for(int iL=0;iL<6;iL++){
00894             if(allRechits[iE][iS][iR][iC][iL].size()){
00895               ++numLayers; 
00896             }
00897           }
00898           if(numLayers>1){
00899             emptyChambers[iE][iS][iR][iC] = false;
00900           }
00901           else{
00902             emptyChambers[iE][iS][iR][iC] = true;
00903           }
00904         }
00905       }
00906     }
00907   }
00908 
00909   //
00910   if (printalot){
00911     printf("  The size of the segment collection is %i\n", int(segments->size()));
00912     //printf("\t...start loop over segments...\n");
00913   }
00914   segmentsPerEvent->Fill(segments->size());
00915   for(CSCSegmentCollection::const_iterator it = segments->begin(); it != segments->end(); it++) {
00916     CSCDetId id  = (CSCDetId)(*it).cscDetId();
00917     StHist[id.endcap()-1][id.station()-1].segmentChi2_ndf->Fill((*it).chi2()/(*it).degreesOfFreedom());
00918     StHist[id.endcap()-1][id.station()-1].hitsInSegment->Fill((*it).nRecHits());
00919     if (printalot){ 
00920       printf("\tendcap/station/ring/chamber: %i %i %i %i\n",
00921              id.endcap(),id.station(),id.ring(),id.chamber());
00922       std::cout<<"\tposition(loc) = "<<(*it).localPosition()<<" error(loc) = "<<(*it).localPositionError()<<std::endl;
00923       std::cout<<"\t chi2/ndf = "<<(*it).chi2()/(*it).degreesOfFreedom()<<" nhits = "<<(*it).nRecHits() <<std::endl;
00924 
00925     }
00926     allSegments[id.endcap()-1][id.station()-1][id.ring()-1][id.chamber()-FirstCh].push_back
00927       (make_pair((*it).localPosition(), (*it).localDirection()));
00928 
00929 
00930     //---- try to get the CSC recHits that contribute to this segment.
00931     //if (printalot) printf("\tGet the recHits for this segment.\t");
00932     std::vector<CSCRecHit2D> theseRecHits = (*it).specificRecHits();
00933     int nRH = (*it).nRecHits();
00934     if (printalot){
00935       printf("\tGet the recHits for this segment.\t");
00936       printf("    nRH = %i\n",nRH);
00937     }
00938     //---- Find which of the rechits in the chamber is in the segment
00939     int layerRH = 0;
00940     for ( vector<CSCRecHit2D>::const_iterator iRH = theseRecHits.begin(); iRH != theseRecHits.end(); iRH++) {
00941       ++layerRH;
00942       CSCDetId idRH = (CSCDetId)(*iRH).cscDetId();
00943       if(printalot){ 
00944         printf("\t%i RH\tendcap/station/ring/chamber/layer: %i/%i/%i/%i/%i\n",
00945                layerRH,idRH.endcap(),idRH.station(),idRH.ring(),idRH.chamber(),idRH.layer());
00946       }
00947       for(size_t jRH = 0; 
00948           jRH<allRechits[idRH.endcap()-1][idRH.station()-1][idRH.ring()-1][idRH.chamber()-FirstCh][idRH.layer()-1].size();
00949           ++jRH){
00950         allRechits[idRH.endcap()-1][idRH.station()-1][idRH.ring()-1][idRH.chamber()-FirstCh][idRH.layer()-1][jRH].first;
00951         float xDiff = iRH->localPosition().x() - 
00952           allRechits[idRH.endcap()-1][idRH.station()-1][idRH.ring()-1][idRH.chamber()-FirstCh][idRH.layer()-1][jRH].first.x();
00953         float yDiff = iRH->localPosition().y() - 
00954           allRechits[idRH.endcap()-1][idRH.station()-1][idRH.ring()-1][idRH.chamber()-FirstCh][idRH.layer()-1][jRH].first.y();
00955         if(fabs(xDiff)<0.0001 && fabs(yDiff)<0.0001){
00956           std::pair <LocalPoint, bool> 
00957             recHitPos(allRechits[idRH.endcap()-1][idRH.station()-1][idRH.ring()-1][idRH.chamber()-FirstCh][idRH.layer()-1][jRH].first, true);
00958           allRechits[idRH.endcap()-1][idRH.station()-1][idRH.ring()-1][idRH.chamber()-FirstCh][idRH.layer()-1][jRH] = recHitPos;
00959           if(printalot){ 
00960             std::cout<<" number of the rechit (from zero) in the segment = "<< jRH<<std::endl;
00961           }
00962         }
00963       }
00964     }
00965   }
00966 }
00967 //
00968 
00969 void CSCEfficiency::ringCandidates(int station, float feta, std::map <std::string, bool> & chamberTypes){
00970   // yeah, hardcoded again...
00971   switch (station){
00972   case 1:
00973     if(feta>0.85 && feta<1.18){//ME13
00974       chamberTypes["ME13"] = true;
00975     }      
00976     if(feta>1.18 && feta<1.7){//ME12
00977       chamberTypes["ME12"] = true;
00978     }
00979     if(feta>1.5 && feta<2.45){//ME11
00980       chamberTypes["ME11"] = true;
00981     }
00982     break;
00983   case 2:
00984     if(feta>0.95 && feta<1.6){//ME22
00985       chamberTypes["ME22"] = true;
00986       
00987     }  
00988     if(feta>1.55 && feta<2.45){//ME21
00989       chamberTypes["ME21"] = true;
00990     }
00991     break;
00992   case 3:
00993     if(feta>1.08 && feta<1.72){//ME32
00994       chamberTypes["ME32"] = true;
00995       
00996     }  
00997     if(feta>1.69 && feta<2.45){//ME31
00998       chamberTypes["ME31"] = true;
00999     }
01000     break;
01001   case 4:
01002     if(feta>1.78 && feta<2.45){//ME41
01003       chamberTypes["ME41"] = true;
01004     }
01005     break;
01006   default:
01007     break;
01008   }
01009 }
01010 //
01011 void CSCEfficiency::chamberCandidates(int station, int ring, float phi, std::vector <int> &coupleOfChambers){
01012   coupleOfChambers.clear();
01013   // -pi< phi<+pi
01014   float phi_zero = 0.;// check! the phi at the "edge" of Ch 1
01015   float phi_const = 2.*M_PI/36.;
01016   int last_chamber = 36;
01017   int first_chamber = 1;
01018   if(1 != station && 1==ring){ // 18 chambers in the ring
01019     phi_const*=2;
01020     last_chamber /= 2;
01021   }
01022   if(phi<0.){
01023     if (printalot) std::cout<<" info: negative phi = "<<phi<<std::endl;
01024     phi += 2*M_PI;
01025   }
01026   float chamber_float = (phi - phi_zero)/phi_const;
01027   int chamber_int = int(chamber_float);
01028   if (chamber_float - float(chamber_int) -0.5 <0.){
01029     if(0!=chamber_int ){
01030       coupleOfChambers.push_back(chamber_int);
01031     }
01032     else{
01033       coupleOfChambers.push_back(last_chamber);
01034     }
01035     coupleOfChambers.push_back(chamber_int+1); 
01036     
01037   }
01038   else{
01039     coupleOfChambers.push_back(chamber_int+1);
01040     if(last_chamber!=chamber_int+1){
01041       coupleOfChambers.push_back(chamber_int+2);
01042     } 
01043     else{
01044       coupleOfChambers.push_back(first_chamber);
01045     }
01046   }
01047   if (printalot) std::cout<<" phi = "<<phi<<" phi_zero = "<<phi_zero<<" phi_const = "<<phi_const<<
01048     " candidate chambers: first ch = "<<coupleOfChambers[0]<<" second ch = "<<coupleOfChambers[1]<<std::endl;
01049 }
01050 
01051 //
01052 bool CSCEfficiency::efficienciesPerChamber(CSCDetId & id, const CSCChamber* cscChamber, FreeTrajectoryState &ftsChamber){
01053   int ec, st, rg, ch, secondRing;
01054   returnTypes(id, ec, st, rg, ch, secondRing);
01055 
01056   LocalVector localDir = cscChamber->toLocal(ftsChamber.momentum());
01057   if(printalot){
01058     std::cout<<" global dir = "<<ftsChamber.momentum()<<std::endl;
01059     std::cout<<" local dir = "<<localDir<<std::endl;
01060     std::cout<<" local theta = "<<localDir.theta()<<std::endl;
01061   }
01062   float dxdz = localDir.x()/localDir.z();
01063   float dydz = localDir.y()/localDir.z();
01064   if(2==st || 3==st){
01065     if(printalot){
01066       std::cout<<"st 3 or 4 ... flip dy/dz"<<std::endl;
01067     }
01068     dydz = - dydz;
01069   }
01070   if(printalot){
01071     std::cout<<"dy/dz = "<<dydz<<std::endl;
01072   }
01073   // Apply angle cut
01074   bool out = true;
01075   if(applyIPangleCuts){
01076     if(dydz>local_DY_DZ_Max || dydz<local_DY_DZ_Min || fabs(dxdz)>local_DX_DZ_Max){ 
01077       out = false;
01078     }
01079   }
01080 
01081   // Segments
01082   bool firstCondition = allSegments[ec][st][rg][ch].size() ? true : false;
01083   bool secondCondition = false;
01084   //---- ME1 is special as usual - ME1a and ME1b are actually one chamber
01085   if(secondRing>-1){
01086     secondCondition = allSegments[ec][st][secondRing][ch].size() ? true : false;
01087   }
01088   if(firstCondition || secondCondition){
01089     if(out){
01090       ChHist[ec][st][rg][ch].digiAppearanceCount->Fill(1);
01091     }  
01092   }
01093   else{
01094     if(out){
01095       ChHist[ec][st][rg][ch].digiAppearanceCount->Fill(0);
01096     } 
01097   }
01098 
01099   if(useDigis){
01100     // ALCTs
01101     firstCondition = allALCT[ec][st][rg][ch];
01102     secondCondition = false;
01103     if(secondRing>-1){
01104       secondCondition = allALCT[ec][st][secondRing][ch];
01105     }
01106     if(firstCondition || secondCondition){
01107       if(out){
01108         ChHist[ec][st][rg][ch].digiAppearanceCount->Fill(3);
01109       }
01110       // always apply partial angle cuts for this kind of histos 
01111       if(fabs(dxdz)<local_DX_DZ_Max){
01112         StHist[ec][st].EfficientALCT_momTheta->Fill(ftsChamber.momentum().theta());
01113         ChHist[ec][st][rg][ch].EfficientALCT_dydz->Fill(dydz);
01114       }
01115     }
01116     else{
01117       if(out){
01118         ChHist[ec][st][rg][ch].digiAppearanceCount->Fill(2);
01119       }
01120       if(fabs(dxdz)<local_DX_DZ_Max){ 
01121         StHist[ec][st].InefficientALCT_momTheta->Fill(ftsChamber.momentum().theta());
01122         ChHist[ec][st][rg][ch].InefficientALCT_dydz->Fill(dydz);
01123       }
01124       if(printalot){
01125         std::cout<<" missing ALCT (dy/dz = "<<dydz<<")";
01126         printf("\t\tendcap/station/ring/chamber: %i/%i/%i/%i\n",ec+1,st+1,rg+1,ch+1);
01127       }
01128     }
01129     
01130     // CLCTs
01131     firstCondition = allCLCT[ec][st][rg][ch];
01132     secondCondition = false;
01133     if(secondRing>-1){
01134       secondCondition = allCLCT[ec][st][secondRing][ch];
01135     }
01136     if(firstCondition || secondCondition){
01137       if(out){
01138         ChHist[ec][st][rg][ch].digiAppearanceCount->Fill(5);
01139       }
01140       if(dydz<local_DY_DZ_Max && dydz>local_DY_DZ_Min){
01141         StHist[ec][st].EfficientCLCT_momPhi->Fill(ftsChamber.momentum().phi() );// - phi chamber...
01142         ChHist[ec][st][rg][ch].EfficientCLCT_dxdz->Fill(dxdz);
01143       }
01144     }
01145     else{
01146       if(out){
01147         ChHist[ec][st][rg][ch].digiAppearanceCount->Fill(4);
01148       }
01149       if(dydz<local_DY_DZ_Max && dydz>local_DY_DZ_Min){ 
01150         StHist[ec][st].InefficientCLCT_momPhi->Fill(ftsChamber.momentum().phi());// - phi chamber...
01151         ChHist[ec][st][rg][ch].InefficientCLCT_dxdz->Fill(dxdz);
01152       }
01153       if(printalot){
01154         std::cout<<" missing CLCT  (dx/dz = "<<dxdz<<")";
01155         printf("\t\tendcap/station/ring/chamber: %i/%i/%i/%i\n",ec+1,st+1,rg+1,ch+1);
01156       }
01157     }
01158     if(out){
01159       // CorrLCTs
01160       firstCondition = allCorrLCT[ec][st][rg][ch];
01161       secondCondition = false;
01162       if(secondRing>-1){
01163         secondCondition = allCorrLCT[ec][st][secondRing][ch];
01164       }
01165       if(firstCondition || secondCondition){
01166         ChHist[ec][st][rg][ch].digiAppearanceCount->Fill(7);
01167       }
01168       else{
01169         ChHist[ec][st][rg][ch].digiAppearanceCount->Fill(6);
01170       }
01171     }
01172   }
01173   return out;
01174 }
01175 
01176 //
01177 bool CSCEfficiency::stripWire_Efficiencies(CSCDetId & id, FreeTrajectoryState &ftsChamber){
01178   int ec, st, rg, ch, secondRing;
01179   returnTypes(id, ec, st, rg, ch, secondRing);
01180 
01181   bool firstCondition, secondCondition;
01182   int missingLayers_s = 0;
01183   int missingLayers_wg = 0;
01184   for(int iLayer=0;iLayer<6;iLayer++){
01185  //----Strips
01186     if(printalot){ 
01187       printf("\t%i swEff: \tendcap/station/ring/chamber/layer: %i/%i/%i/%i/%i\n",
01188              iLayer + 1,id.endcap(),id.station(),id.ring(),id.chamber(),iLayer+1);
01189       std::cout<<" size S = "<<allStrips[id.endcap()-1][id.station()-1][id.ring()-1][id.chamber()-FirstCh][iLayer].size()<<
01190         "size W = "<<allWG[id.endcap()-1][id.station()-1][id.ring()-1][id.chamber()-FirstCh][iLayer].size()<<std::endl;
01191 
01192     }
01193     firstCondition = allStrips[ec][st][rg][ch][iLayer].size() ? true : false;
01194     //allSegments[ec][st][rg][ch].size() ? true : false;
01195     secondCondition = false;
01196     if(secondRing>-1){
01197       secondCondition = allStrips[ec][st][secondRing][ch][iLayer].size() ? true : false;
01198     }
01199     if(firstCondition || secondCondition){
01200       ChHist[ec][st][rg][ch].EfficientStrips->Fill(iLayer+1);
01201     }
01202     else{
01203       if(printalot){ 
01204         std::cout<<"missing strips ";
01205         printf("\t\tendcap/station/ring/chamber/layer: %i/%i/%i/%i/%i\n",id.endcap(),id.station(),id.ring(),id.chamber(),iLayer+1);
01206       }
01207     }
01208     // Wires
01209     firstCondition = allWG[ec][st][rg][ch][iLayer].size() ? true : false;
01210     secondCondition = false;
01211     if(secondRing>-1){
01212       secondCondition = allWG[ec][st][secondRing][ch][iLayer].size() ? true : false;
01213     }
01214     if(firstCondition || secondCondition){
01215       ChHist[ec][st][rg][ch].EfficientWireGroups->Fill(iLayer+1);
01216     }
01217     else{
01218       if(printalot){ 
01219         std::cout<<"missing wires ";
01220         printf("\t\tendcap/station/ring/chamber/layer: %i/%i/%i/%i/%i\n",id.endcap(),id.station(),id.ring(),id.chamber(),iLayer+1);
01221       }
01222     }
01223   }
01224   // Normalization
01225   if(6!=missingLayers_s){
01226     ChHist[ec][st][rg][ch].EfficientStrips->Fill(8);
01227   }
01228   if(6!=missingLayers_wg){
01229     ChHist[ec][st][rg][ch].EfficientWireGroups->Fill(8);
01230   }
01231   ChHist[ec][st][rg][ch].EfficientStrips->Fill(9);
01232   ChHist[ec][st][rg][ch].EfficientWireGroups->Fill(9);
01233 //
01234   ChHist[ec][st][rg][ch].StripWiresCorrelations->Fill(1);
01235   if(missingLayers_s!=missingLayers_wg){
01236     ChHist[ec][st][rg][ch].StripWiresCorrelations->Fill(2);
01237     if(6==missingLayers_wg){
01238       ChHist[ec][st][rg][ch].StripWiresCorrelations->Fill(3);
01239       ChHist[ec][st][rg][ch].NoWires_momTheta->Fill(ftsChamber.momentum().theta());
01240     }
01241     if(6==missingLayers_s){
01242       ChHist[ec][st][rg][ch].StripWiresCorrelations->Fill(4);
01243       ChHist[ec][st][rg][ch].NoStrips_momPhi->Fill(ftsChamber.momentum().theta());
01244     }
01245   }
01246   else if(6==missingLayers_s){
01247     ChHist[ec][st][rg][ch].StripWiresCorrelations->Fill(5);
01248   }
01249 
01250   return true;  
01251 }
01252 //
01253 bool CSCEfficiency::recSimHitEfficiency(CSCDetId & id, FreeTrajectoryState &ftsChamber){
01254   int ec, st, rg, ch, secondRing;
01255   returnTypes(id, ec, st, rg, ch, secondRing);
01256   bool firstCondition, secondCondition;
01257   for(int iLayer=0; iLayer<6;iLayer++){
01258     firstCondition = allSimhits[ec][st][rg][ch][iLayer].size() ? true : false;
01259     secondCondition = false;
01260     int thisRing = rg;
01261     if(secondRing>-1){
01262       secondCondition = allSimhits[ec][st][secondRing][ch][iLayer].size() ? true : false;
01263       if(secondCondition){
01264         thisRing = secondRing;
01265       }
01266     }
01267     if(firstCondition || secondCondition){
01268       for(size_t iSH=0;
01269           iSH<allSimhits[ec][st][thisRing][ch][iLayer].size();
01270           ++iSH){
01271         if(13 ==
01272            fabs(allSimhits[ec][st][thisRing][ch][iLayer][iSH].second)){
01273           ChHist[ec][st][rg][ch].SimSimhits->Fill(iLayer+1);
01274           if(allRechits[ec][st][thisRing][ch][iLayer].size()){
01275             ChHist[ec][st][rg][ch].SimRechits->Fill(iLayer+1);
01276           }
01277           break;
01278         }
01279       }
01280       //---- Next is not too usefull... 
01281       /*
01282       for(unsigned int iSimHits=0;
01283           iSimHits<allSimhits[id.endcap()-1][id.station()-1][id.ring()-1][id.chamber()-FirstCh][iLayer].size();
01284           iSimHits++){
01285         ChHist[ec][st][rg][id.chamber()-FirstCh].SimSimhits_each->Fill(iLayer+1);
01286       }
01287       for(unsigned int iRecHits=0;
01288           iRecHits<allRechits[id.endcap()-1][id.station()-1][id.ring()-1][id.chamber()-FirstCh][iLayer].size();
01289           iRecHits++){
01290         ChHist[ec][st][rg][id.chamber()-FirstCh].SimRechits_each->Fill(iLayer+1);
01291       }
01292       */
01293       //
01294     }   
01295   }
01296   return true;
01297 }
01298 
01299 //
01300 bool CSCEfficiency::recHitSegment_Efficiencies(CSCDetId & id, const CSCChamber* cscChamber, FreeTrajectoryState &ftsChamber){
01301   int ec, st, rg, ch, secondRing;
01302   returnTypes(id, ec, st, rg, ch, secondRing);
01303   bool firstCondition, secondCondition;
01304 
01305   std::vector <bool> missingLayers_rh(6);
01306   std::vector <int> usedInSegment(6);
01307   // Rechits
01308   if(printalot) std::cout<<"RecHits eff"<<std::endl;
01309   for(int iLayer=0;iLayer<6;++iLayer){
01310     firstCondition = allRechits[ec][st][rg][ch][iLayer].size() ? true : false;
01311     secondCondition = false;
01312     int thisRing = rg;
01313     if(secondRing>-1){
01314       secondCondition = allRechits[ec][st][secondRing][ch][iLayer].size() ? true : false;
01315       if(secondCondition){
01316         thisRing = secondRing;
01317       }
01318     }
01319     if(firstCondition || secondCondition){
01320       ChHist[ec][st][rg][ch].EfficientRechits_good->Fill(iLayer+1);
01321       for(size_t iR=0;
01322           iR<allRechits[ec][st][thisRing][ch][iLayer].size();
01323           ++iR){
01324         if(allRechits[ec][st][thisRing][ch][iLayer][iR].second){
01325           usedInSegment[iLayer] = 1;
01326           break;
01327         }
01328         else{
01329           usedInSegment[iLayer] = -1;
01330         }
01331       }
01332     }
01333     else{
01334       missingLayers_rh[iLayer] = true;
01335       if(printalot){
01336         std::cout<<"missing rechits ";
01337         printf("\t\tendcap/station/ring/chamber/layer: %i/%i/%i/%i/%i\n",id.endcap(),id.station(),id.ring(),id.chamber(),iLayer+1);
01338       }
01339     }
01340   }
01341   GlobalVector globalDir;
01342   GlobalPoint globalPos;
01343  // Segments
01344   firstCondition = allSegments[ec][st][rg][ch].size() ? true : false;
01345   secondCondition = false;
01346   int secondSize = 0;
01347   int thisRing = rg;
01348   if(secondRing>-1){
01349     secondCondition = allSegments[ec][st][secondRing][ch].size() ? true : false;
01350     secondSize = allSegments[ec][st][secondRing][ch].size();
01351     if(secondCondition){
01352       thisRing = secondRing;
01353     }
01354   }
01355   if(firstCondition || secondCondition){
01356     if (printalot) std::cout<<"segments - start ec = "<<ec<<" st = "<<st<<" rg = "<<rg<<" ch = "<<ch<<std::endl;
01357     StHist[ec][st].EfficientSegments_XY->Fill(ftsChamber.position().x(),ftsChamber.position().y());
01358     if(1==allSegments[ec][st][rg][ch].size() + secondSize){
01359       globalDir = cscChamber->toGlobal(allSegments[ec][st][thisRing][ch][0].second);
01360       globalPos = cscChamber->toGlobal(allSegments[ec][st][thisRing][ch][0].first);
01361       StHist[ec][st].EfficientSegments_eta->Fill(fabs(ftsChamber.position().eta()));
01362       double residual = sqrt(pow(ftsChamber.position().x() - globalPos.x(),2)+
01363                              pow(ftsChamber.position().y() - globalPos.y(),2)+
01364                              pow(ftsChamber.position().z() - globalPos.z(),2));
01365       if (printalot) std::cout<<" fts.position() = "<<ftsChamber.position()<<" segPos = "<<globalPos<<" res = "<<residual<< std::endl;
01366       StHist[ec][st].ResidualSegments->Fill(residual);
01367     }
01368     for(int iLayer=0;iLayer<6;++iLayer){
01369       if(printalot) std::cout<<" iLayer = "<<iLayer<<" usedInSegment = "<<usedInSegment[iLayer]<<std::endl;
01370       if(0!=usedInSegment[iLayer]){
01371         if(-1==usedInSegment[iLayer]){
01372           ChHist[ec][st][rg][ch].InefficientSingleHits->Fill(iLayer+1);
01373         }
01374         ChHist[ec][st][rg][ch].AllSingleHits->Fill(iLayer+1);
01375       }
01376       firstCondition = allRechits[ec][st][rg][ch][iLayer].size() ? true : false;
01377       secondCondition = false;
01378       if(secondRing>-1){
01379         secondCondition = allRechits[ec][st][secondRing][ch][iLayer].size() ? true : false;
01380       }
01381       float stripAngle = 99999.;
01382       std::vector<float> posXY(2);
01383       bool oneSegment = false;
01384       if(1==allSegments[ec][st][rg][ch].size() + secondSize){
01385         oneSegment = true;
01386         const BoundPlane bp = cscChamber->layer(iLayer+1)->surface();
01387         linearExtrapolation(globalPos,globalDir, bp.position().z(), posXY);
01388         GlobalPoint gp_extrapol( posXY.at(0), posXY.at(1),bp.position().z());
01389         const LocalPoint lp_extrapol = cscChamber->layer(iLayer+1)->toLocal(gp_extrapol);
01390         posXY.at(0) = lp_extrapol.x();
01391         posXY.at(1) = lp_extrapol.y();
01392         int nearestStrip = cscChamber->layer(iLayer+1)->geometry()->nearestStrip(lp_extrapol);
01393         stripAngle = cscChamber->layer(iLayer+1)->geometry()->stripAngle(nearestStrip) - M_PI/2. ;
01394       }
01395       if(firstCondition || secondCondition){
01396         ChHist[ec][st][rg][ch].EfficientRechits_inSegment->Fill(iLayer+1);
01397         if(oneSegment){
01398           ChHist[ec][st][rg][ch].Y_EfficientRecHits_inSegment[iLayer]->Fill(posXY.at(1));
01399           ChHist[ec][st][rg][ch].Phi_EfficientRecHits_inSegment[iLayer]->Fill(stripAngle);
01400         }
01401       }
01402       else{
01403         if(oneSegment){
01404           ChHist[ec][st][rg][ch].Y_InefficientRecHits_inSegment[iLayer]->Fill(posXY.at(1));
01405           ChHist[ec][st][rg][ch].Phi_InefficientRecHits_inSegment[iLayer]->Fill(stripAngle);
01406         }
01407       }
01408     }
01409   }
01410   else{
01411     StHist[ec][st].InefficientSegments_XY->Fill(ftsChamber.position().x(),ftsChamber.position().y());
01412     if(printalot){
01413       std::cout<<"missing segment "<<std::endl;
01414       printf("\t\tendcap/station/ring/chamber: %i/%i/%i/%i\n",id.endcap(),id.station(),id.ring(),id.chamber());
01415       std::cout<<" fts.position() = "<<ftsChamber.position()<<std::endl;
01416     }
01417   }
01418   // Normalization
01419   ChHist[ec][st][rg][ch].EfficientRechits_good->Fill(8);
01420   if(allSegments[ec][st][rg][ch].size()+secondSize<2){
01421     StHist[ec][st].AllSegments_eta->Fill(fabs(ftsChamber.position().eta()));
01422   }
01423   ChHist[ec][st][rg][id.chamber()-FirstCh].EfficientRechits_inSegment->Fill(9);
01424 
01425   return true;
01426 }
01427 //
01428 void CSCEfficiency::returnTypes(CSCDetId & id, int &ec, int &st, int &rg, int &ch, int &secondRing){
01429   ec = id.endcap()-1;
01430   st = id.station()-1;
01431   rg = id.ring()-1;
01432   secondRing = -1;
01433   if(1==id.station() && (4==id.ring() || 1==id.ring()) ){
01434     rg = 0;
01435     secondRing = 3;
01436   }
01437   ch = id.chamber()-FirstCh;
01438 }
01439 
01440 //
01441 void CSCEfficiency::getFromFTS( const FreeTrajectoryState& fts,
01442                                 CLHEP::Hep3Vector& p3, CLHEP::Hep3Vector& r3,
01443                                 int& charge, AlgebraicSymMatrix66& cov){
01444   
01445   GlobalVector p3GV = fts.momentum();
01446   GlobalPoint r3GP = fts.position();
01447 
01448   p3.set(p3GV.x(), p3GV.y(), p3GV.z());
01449   r3.set(r3GP.x(), r3GP.y(), r3GP.z());
01450 
01451   charge = fts.charge();
01452   cov = fts.hasError() ? fts.cartesianError().matrix() : AlgebraicSymMatrix66();
01453 
01454 }
01455 
01456 FreeTrajectoryState CSCEfficiency::getFromCLHEP(const CLHEP::Hep3Vector& p3, const CLHEP::Hep3Vector& r3,
01457                                                        int charge, const AlgebraicSymMatrix66& cov,
01458                                                        const MagneticField* field){
01459 
01460   GlobalVector p3GV(p3.x(), p3.y(), p3.z());
01461   GlobalPoint r3GP(r3.x(), r3.y(), r3.z());
01462   GlobalTrajectoryParameters tPars(r3GP, p3GV, charge, field);
01463 
01464   CartesianTrajectoryError tCov(cov);
01465 
01466   return cov.kRows == 6 ? FreeTrajectoryState(tPars, tCov) : FreeTrajectoryState(tPars) ;
01467 }
01468 
01469 void CSCEfficiency::linearExtrapolation(GlobalPoint initialPosition ,GlobalVector initialDirection, 
01470                                         float zSurface, std::vector <float> &posZY){
01471   double paramLine =  lineParameter(initialPosition.z(), zSurface, initialDirection.z());
01472   double xPosition = extrapolate1D(initialPosition.x(), initialDirection.x(),paramLine);
01473   double yPosition = extrapolate1D(initialPosition.y(), initialDirection.y(),paramLine);
01474   posZY.clear();
01475   posZY.push_back(xPosition);
01476   posZY.push_back(yPosition);
01477 }
01478 //
01479 double CSCEfficiency::extrapolate1D(double initPosition, double initDirection, double parameterOfTheLine){
01480   double extrapolatedPosition = initPosition + initDirection*parameterOfTheLine;
01481   return extrapolatedPosition;
01482 }
01483 //
01484 double CSCEfficiency::lineParameter(double initZPosition, double destZPosition, double initZDirection){
01485   double paramLine = (destZPosition-initZPosition)/initZDirection;
01486   return paramLine;
01487 }
01488 //
01489 void CSCEfficiency::chooseDirection(CLHEP::Hep3Vector & innerPosition, CLHEP::Hep3Vector & outerPosition){
01490 
01491   //---- Be careful with trigger conditions too
01492   if(!isIPdata){
01493     float dy = outerPosition.y() - innerPosition.y();
01494     float dz = outerPosition.z() - innerPosition.z();
01495     if(isBeamdata){
01496       if(dz>0){
01497         alongZ = true;
01498       }
01499       else{
01500         alongZ = false;
01501       }
01502     }
01503     else{//cosmics
01504       if(dy/dz>0){
01505         alongZ = false;
01506       }
01507       else{
01508         alongZ = true;
01509       }
01510     }
01511   }
01512 }
01513 //
01514 const Propagator* CSCEfficiency::propagator(std::string propagatorName) const {
01515   return &*theService->propagator(propagatorName);
01516 }
01517   
01518 //
01519 TrajectoryStateOnSurface CSCEfficiency::propagate(FreeTrajectoryState & ftsStart, const BoundPlane &bpDest){
01520   TrajectoryStateOnSurface tSOSDest;
01521   std::string propagatorName;
01522 /*
01523 // it would work if cosmic muons had properly assigned direction...
01524   bool dzPositive = bpDest.position().z() - ftsStart.position().z() > 0 ? true : false;
01525  //---- Be careful with trigger conditions too
01526   if(!isIPdata){
01527     bool rightDirection = !(alongZ^dzPositive);
01528     if(rightDirection){
01529       if(printalot) std::cout<<" propagate along momentum"<<std::endl;
01530       propagatorName = "SteppingHelixPropagatorAlong";
01531     }
01532     else{
01533       if(printalot) std::cout<<" propagate opposite momentum"<<std::endl;
01534       propagatorName = "SteppingHelixPropagatorOpposite";
01535     }
01536   }
01537   else{
01538     if(printalot) std::cout<<" propagate any (momentum)"<<std::endl;
01539     propagatorName = "SteppingHelixPropagatorAny";
01540   }
01541 */
01542   propagatorName = "SteppingHelixPropagatorAny";
01543   tSOSDest = propagator(propagatorName)->propagate(ftsStart, bpDest);   
01544   return tSOSDest;
01545 }
01546 //
01547 bool CSCEfficiency::applyTrigger(edm::Handle<edm::TriggerResults> &hltR,
01548                                  const edm::TriggerNames & triggerNames){
01549   bool triggerPassed = true;
01550   std::vector<std::string>  hlNames=triggerNames.triggerNames();
01551   pointToTriggers.clear();
01552   for(size_t imyT = 0;imyT<myTriggers.size();++imyT){
01553     for (size_t iT=0; iT<hlNames.size(); ++iT) {
01554       //std::cout<<" iT = "<<iT<<" hlNames[iT] = "<<hlNames[iT]<<
01555       //" : wasrun = "<<hltR->wasrun(iT)<<" accept = "<<
01556       //         hltR->accept(iT)<<" !error = "<< 
01557       //        !hltR->error(iT)<<std::endl; 
01558       if(!imyT){       
01559         if(hltR->wasrun(iT) &&
01560            hltR->accept(iT) &&
01561            !hltR->error(iT) ){
01562            TriggersFired->Fill(iT); 
01563         }
01564       }
01565       if(hlNames[iT]==myTriggers[imyT]){
01566         pointToTriggers.push_back(iT);
01567         if(imyT){
01568           break;
01569         }
01570       }
01571     }
01572   }
01573   if(pointToTriggers.size()!=myTriggers.size()){
01574     pointToTriggers.clear();
01575     if(printalot){
01576       std::cout<<" Not all trigger names found - all trigger specifications will be ignored. Check your cfg file!"<<std::endl;
01577     }
01578   }
01579   else{
01580     if(pointToTriggers.size()){
01581       if(printalot){
01582         std::cout<<"The following triggers will be required in the event: "<<std::endl;
01583         for(size_t imyT =0; imyT <pointToTriggers.size();++imyT){
01584           std::cout<<"  "<<hlNames[pointToTriggers[imyT]];
01585         }
01586         std::cout<<std::endl;
01587         std::cout<<" in condition (AND/OR) : "<<!andOr<<"/"<<andOr<<std::endl;
01588       }
01589     }
01590   }
01591 
01592   if (hltR.isValid()) {
01593     if(!pointToTriggers.size()){
01594       if(printalot){
01595         std::cout<<" No triggers specified in the configuration or all ignored - no trigger information will be considered"<<std::endl;
01596       }
01597     }
01598     for(size_t imyT =0; imyT <pointToTriggers.size();++imyT){
01599       if(hltR->wasrun(pointToTriggers[imyT]) && 
01600          hltR->accept(pointToTriggers[imyT]) && 
01601          !hltR->error(pointToTriggers[imyT]) ){
01602         triggerPassed = true;
01603         if(andOr){
01604           break;
01605         }
01606       }
01607       else{
01608         triggerPassed = false;
01609         if(!andOr){
01610           triggerPassed = false;
01611           break;
01612         }
01613       }
01614     }
01615   }
01616   else{
01617     if(printalot){
01618       std::cout<<" TriggerResults handle returns invalid state?! No trigger information will be considered"<<std::endl;
01619     }
01620   }
01621   if(printalot){
01622     std::cout<<" Trigger passed: "<<triggerPassed<<std::endl;
01623   }
01624   return triggerPassed;
01625 }
01626 //
01627 
01628 // Constructor
01629 CSCEfficiency::CSCEfficiency(const ParameterSet& pset){
01630 
01631   // const float Xmin = -70;
01632   //const float Xmax = 70;
01633   //const int nXbins = int(4.*(Xmax - Xmin));
01634   const float Ymin = -165;
01635   const float Ymax = 165;
01636   const int nYbins = int((Ymax - Ymin)/2);
01637   const float Layer_min = -0.5;
01638   const float Layer_max = 9.5;
01639   const int nLayer_bins = int(Layer_max - Layer_min);
01640   //
01641 
01642   //---- Get the input parameters
01643   printout_NEvents  = pset.getUntrackedParameter<unsigned int>("printout_NEvents",0);
01644   rootFileName     = pset.getUntrackedParameter<string>("rootFileName","cscHists.root");
01645 
01646   isData  = pset.getUntrackedParameter<bool>("runOnData",true);// 
01647   isIPdata  = pset.getUntrackedParameter<bool>("IPdata",false);// 
01648   isBeamdata  = pset.getUntrackedParameter<bool>("Beamdata",false);//
01649   getAbsoluteEfficiency  = pset.getUntrackedParameter<bool>("getAbsoluteEfficiency",true);//
01650   useDigis = pset.getUntrackedParameter<bool>("useDigis", true);// 
01651   distanceFromDeadZone = pset.getUntrackedParameter<double>("distanceFromDeadZone", 10.);// 
01652   minP = pset.getUntrackedParameter<double>("minP",20.);//
01653   maxP = pset.getUntrackedParameter<double>("maxP",100.);//
01654   maxNormChi2 = pset.getUntrackedParameter<double>("maxNormChi2", 3.);//
01655   minTrackHits = pset.getUntrackedParameter<unsigned int>("minTrackHits",10);//
01656 
01657   applyIPangleCuts = pset.getUntrackedParameter<bool>("applyIPangleCuts", false);// 
01658     local_DY_DZ_Max = pset.getUntrackedParameter<double>("local_DY_DZ_Max",-0.1);//
01659       local_DY_DZ_Min = pset.getUntrackedParameter<double>("local_DY_DZ_Min",-0.8);//
01660         local_DX_DZ_Max = pset.getUntrackedParameter<double>("local_DX_DZ_Max",0.2);//
01661 
01662   alctDigiTag_  = pset.getParameter<edm::InputTag>("alctDigiTag") ;
01663   clctDigiTag_  = pset.getParameter<edm::InputTag>("clctDigiTag") ;
01664   corrlctDigiTag_  = pset.getParameter<edm::InputTag>("corrlctDigiTag") ;
01665   stripDigiTag_  = pset.getParameter<edm::InputTag>("stripDigiTag") ;
01666   wireDigiTag_ = pset.getParameter<edm::InputTag>("wireDigiTag") ;
01667   rechitDigiTag_ = pset.getParameter<edm::InputTag>("rechitDigiTag") ;
01668   segmentDigiTag_ = pset.getParameter<edm::InputTag>("segmentDigiTag") ;
01669   simHitTag     = pset.getParameter<edm::InputTag>("simHitTag");
01670   tracksTag = pset.getParameter< edm::InputTag >("tracksTag");
01671 
01672   ParameterSet serviceParameters = pset.getParameter<ParameterSet>("ServiceParameters");
01673   // maybe use the service for getting magnetic field, propagators, etc. ...
01674   theService        = new MuonServiceProxy(serviceParameters);
01675 
01676   // Trigger
01677   useTrigger =  pset.getUntrackedParameter<bool>("useTrigger", false);
01678   hlTriggerResults_ = pset.getParameter<edm::InputTag> ("HLTriggerResults");
01679   myTriggers = pset.getParameter<std::vector <std::string> >("myTriggers");
01680   andOr =  pset.getUntrackedParameter<bool>("andOr");
01681   pointToTriggers.clear();
01682 
01683 
01684   //---- set counter to zero
01685   nEventsAnalyzed = 0;
01686   //---- set presence of magnetic field
01687   magField = true;
01688   //
01689   std::string Path = "AllChambers/";
01690   std::string FullName;
01691   //---- File with output histograms 
01692   theFile = new TFile(rootFileName.c_str(), "RECREATE");
01693   theFile->cd();
01694   //---- Book histograms for the analysis
01695   char SpecName[50];
01696   
01697   sprintf(SpecName,"DataFlow"); 
01698   DataFlow =  
01699     new TH1F(SpecName,"Data flow;condition number;entries",40,-0.5,39.5);
01700   //
01701   sprintf(SpecName,"TriggersFired"); 
01702   TriggersFired =
01703     new TH1F(SpecName,"Triggers fired;trigger number;entries",140,-0.5,139.5);
01704   //
01705   int Chan = 50;
01706   float minChan = -0.5;
01707   float maxChan = 49.5;
01708   //
01709   sprintf(SpecName,"ALCTPerEvent");    
01710   ALCTPerEvent = new TH1F(SpecName,"ALCTs per event;N digis;entries",Chan,minChan,maxChan);
01711   //
01712   sprintf(SpecName,"CLCTPerEvent");    
01713   CLCTPerEvent = new TH1F(SpecName,"CLCTs per event;N digis;entries",Chan,minChan,maxChan);
01714   //
01715     sprintf(SpecName,"recHitsPerEvent");    
01716   recHitsPerEvent = new TH1F(SpecName,"RecHits per event;N digis;entries",150,-0.5,149.5);
01717   //
01718   sprintf(SpecName,"segmentsPerEvent");    
01719   segmentsPerEvent = new TH1F(SpecName,"segments per event;N digis;entries",Chan,minChan,maxChan);
01720   //
01721   //---- Book groups of histograms (for any chamber)
01722 
01723   map<std::string,bool>::iterator iter;  
01724   for(int ec = 0;ec<2;++ec){
01725     for(int st = 0;st<4;++st){
01726       theFile->cd();
01727       sprintf(SpecName,"Stations__E%d_S%d",ec+1, st+1);
01728       theFile->mkdir(SpecName);
01729       theFile->cd(SpecName);
01730 
01731       //
01732       sprintf(SpecName,"segmentChi2_ndf_St%d",st+1);
01733       StHist[ec][st].segmentChi2_ndf = 
01734         new TH1F(SpecName,"Chi2/ndf of a segment;chi2/ndf;entries",100,0.,20.);
01735       //
01736       sprintf(SpecName,"hitsInSegment_St%d",st+1);
01737       StHist[ec][st].hitsInSegment = 
01738         new TH1F(SpecName,"Number of hits in a segment;nHits;entries",7,-0.5,6.5);
01739       //
01740       Chan = 170;
01741       minChan = 0.85;
01742       maxChan = 2.55;
01743       //
01744       sprintf(SpecName,"AllSegments_eta_St%d",st+1);
01745       StHist[ec][st].AllSegments_eta = 
01746         new TH1F(SpecName,"All segments in eta;eta;entries",Chan,minChan,maxChan);
01747       //
01748       sprintf(SpecName,"EfficientSegments_eta_St%d",st+1);
01749       StHist[ec][st].EfficientSegments_eta = 
01750         new TH1F(SpecName,"Efficient segments in eta;eta;entries",Chan,minChan,maxChan);
01751       //
01752       sprintf(SpecName,"ResidualSegments_St%d",st+1);
01753       StHist[ec][st].ResidualSegments = 
01754         new TH1F(SpecName,"Residual (segments);residual,cm;entries",75,0.,15.);
01755       //
01756       Chan = 200;
01757       minChan = -800.;
01758       maxChan = 800.;
01759       int Chan2 = 200;
01760       float minChan2 = -800.;
01761       float maxChan2 = 800.;
01762       
01763       sprintf(SpecName,"EfficientSegments_XY_St%d",st+1);
01764       StHist[ec][st].EfficientSegments_XY = new TH2F(SpecName,"Efficient segments in XY;X;Y",
01765                                                      Chan,minChan,maxChan,Chan2,minChan2,maxChan2);
01766       sprintf(SpecName,"InefficientSegments_XY_St%d",st+1);
01767       StHist[ec][st].InefficientSegments_XY = new TH2F(SpecName,"Inefficient segments in XY;X;Y",
01768                                                        Chan,minChan,maxChan,Chan2,minChan2,maxChan2);
01769       //
01770       Chan = 80;
01771       minChan = 0;
01772       maxChan = 3.2;
01773       sprintf(SpecName,"EfficientALCT_momTheta_St%d",st+1);
01774       StHist[ec][st].EfficientALCT_momTheta = new TH1F(SpecName,"Efficient ALCT in theta (momentum);theta, rad;entries",
01775                                                        Chan,minChan,maxChan);
01776       //
01777       sprintf(SpecName,"InefficientALCT_momTheta_St%d",st+1);
01778       StHist[ec][st].InefficientALCT_momTheta = new TH1F(SpecName,"Inefficient ALCT in theta (momentum);theta, rad;entries",
01779                                                          Chan,minChan,maxChan);
01780       //
01781       Chan = 160;
01782       minChan = -3.2;
01783       maxChan = 3.2;
01784       sprintf(SpecName,"EfficientCLCT_momPhi_St%d",st+1);
01785       StHist[ec][st].EfficientCLCT_momPhi = new TH1F(SpecName,"Efficient CLCT in phi (momentum);phi, rad;entries",
01786                                                      Chan,minChan,maxChan);
01787       //
01788       sprintf(SpecName,"InefficientCLCT_momPhi_St%d",st+1);
01789       StHist[ec][st].InefficientCLCT_momPhi = new TH1F(SpecName,"Inefficient CLCT in phi (momentum);phi, rad;entries",
01790                                                        Chan,minChan,maxChan);
01791       //
01792       theFile->cd();
01793       for(int rg = 0;rg<3;++rg){
01794         if(0!=st && rg>1){
01795           continue;
01796         }
01797         else if(1==rg && 3==st){
01798           continue;
01799         }
01800         for(int iChamber=FirstCh;iChamber<FirstCh+NumCh;iChamber++){
01801           if(0!=st && 0==rg && iChamber >18){
01802             continue;
01803           }
01804           theFile->cd();
01805           sprintf(SpecName,"Chambers__E%d_S%d_R%d_Chamber_%d",ec+1, st+1, rg+1,iChamber);
01806           theFile->mkdir(SpecName);
01807           theFile->cd(SpecName);
01808           //
01809 
01810           sprintf(SpecName,"EfficientRechits_inSegment_Ch%d",iChamber);
01811           ChHist[ec][st][rg][iChamber-FirstCh].EfficientRechits_inSegment = 
01812             new TH1F(SpecName,"Existing RecHit given a segment;layers (1-6);entries",nLayer_bins,Layer_min,Layer_max);
01813           //
01814           sprintf(SpecName,"InefficientSingleHits_Ch%d",iChamber);
01815           ChHist[ec][st][rg][iChamber-FirstCh].InefficientSingleHits = 
01816             new TH1F(SpecName,"Single RecHits not in the segment;layers (1-6);entries ",nLayer_bins,Layer_min,Layer_max);
01817           //
01818           sprintf(SpecName,"AllSingleHits_Ch%d",iChamber);
01819           ChHist[ec][st][rg][iChamber-FirstCh].AllSingleHits = 
01820             new TH1F(SpecName,"Single RecHits given a segment; layers (1-6);entries",nLayer_bins,Layer_min,Layer_max);
01821           //
01822           sprintf(SpecName,"digiAppearanceCount_Ch%d",iChamber);
01823           ChHist[ec][st][rg][iChamber-FirstCh].digiAppearanceCount = 
01824             new TH1F(SpecName,"Digi appearance (no-yes): segment(0,1), ALCT(2,3), CLCT(4,5), CorrLCT(6,7); digi type;entries",
01825                      8,-0.5,7.5);
01826           //
01827           Chan = 100;
01828           minChan = -1.1;
01829           maxChan = 0.9;
01830           sprintf(SpecName,"EfficientALCT_dydz_Ch%d",iChamber);
01831           ChHist[ec][st][rg][iChamber-FirstCh].EfficientALCT_dydz = 
01832             new TH1F(SpecName,"Efficient ALCT; local dy/dz (ME 3 and 4 flipped);entries",
01833                      Chan, minChan, maxChan);
01834           //
01835           sprintf(SpecName,"InefficientALCT_dydz_Ch%d",iChamber);
01836           ChHist[ec][st][rg][iChamber-FirstCh].InefficientALCT_dydz = 
01837             new TH1F(SpecName,"Inefficient ALCT; local dy/dz (ME 3 and 4 flipped);entries",
01838                      Chan, minChan, maxChan);
01839           //
01840           Chan = 100;
01841           minChan = -1.;
01842           maxChan = 1.0;
01843           sprintf(SpecName,"EfficientCLCT_dxdz_Ch%d",iChamber);
01844           ChHist[ec][st][rg][iChamber-FirstCh].EfficientCLCT_dxdz = 
01845             new TH1F(SpecName,"Efficient CLCT; local dxdz;entries",
01846                      Chan, minChan, maxChan);
01847           //
01848           sprintf(SpecName,"InefficientCLCT_dxdz_Ch%d",iChamber);
01849           ChHist[ec][st][rg][iChamber-FirstCh].InefficientCLCT_dxdz = 
01850             new TH1F(SpecName,"Inefficient CLCT; local dxdz;entries",
01851                      Chan, minChan, maxChan);
01852           //
01853           sprintf(SpecName,"EfficientRechits_good_Ch%d",iChamber);
01854           ChHist[ec][st][rg][iChamber-FirstCh].EfficientRechits_good = 
01855             new TH1F(SpecName,"Existing RecHit - sensitive area only;layers (1-6);entries",nLayer_bins,Layer_min,Layer_max);
01856           //
01857           sprintf(SpecName,"EfficientStrips_Ch%d",iChamber);
01858           ChHist[ec][st][rg][iChamber-FirstCh].EfficientStrips = 
01859             new TH1F(SpecName,"Existing strip;layer (1-6); entries",nLayer_bins,Layer_min,Layer_max);
01860           //
01861           sprintf(SpecName,"EfficientWireGroups_Ch%d",iChamber);
01862           ChHist[ec][st][rg][iChamber-FirstCh].EfficientWireGroups = 
01863             new TH1F(SpecName,"Existing WireGroups;layer (1-6); entries ",nLayer_bins,Layer_min,Layer_max);
01864           //
01865           sprintf(SpecName,"StripWiresCorrelations_Ch%d",iChamber);
01866           ChHist[ec][st][rg][iChamber-FirstCh].StripWiresCorrelations = 
01867             new TH1F(SpecName,"StripWire correlations;; entries ",5,0.5,5.5);
01868           //
01869           Chan = 80;
01870           minChan = 0;
01871           maxChan = 3.2;
01872           sprintf(SpecName,"NoWires_momTheta_Ch%d",iChamber);
01873           ChHist[ec][st][rg][iChamber-FirstCh].NoWires_momTheta = 
01874             new TH1F(SpecName,"No wires (all strips present) - in theta (momentum);theta, rad;entries",
01875                      Chan,minChan,maxChan);
01876           //
01877           Chan = 160;
01878           minChan = -3.2;
01879           maxChan = 3.2;
01880           sprintf(SpecName,"NoStrips_momPhi_Ch%d",iChamber);
01881           ChHist[ec][st][rg][iChamber-FirstCh].NoStrips_momPhi = 
01882             new TH1F(SpecName,"No strips (all wires present) - in phi (momentum);phi, rad;entries",
01883                      Chan,minChan,maxChan);
01884           //
01885           for(int iLayer=0; iLayer<6;iLayer++){
01886             sprintf(SpecName,"Y_InefficientRecHits_inSegment_Ch%d_L%d",iChamber,iLayer);
01887             ChHist[ec][st][rg][iChamber-FirstCh].Y_InefficientRecHits_inSegment.push_back
01888               (new TH1F(SpecName,"Missing RecHit/layer in a segment (local system, whole chamber);Y, cm; entries",
01889                         nYbins,Ymin, Ymax));
01890             //
01891             sprintf(SpecName,"Y_EfficientRecHits_inSegment_Ch%d_L%d",iChamber,iLayer);
01892             ChHist[ec][st][rg][iChamber-FirstCh].Y_EfficientRecHits_inSegment.push_back
01893               (new TH1F(SpecName,"Efficient (extrapolated from the segment) RecHit/layer in a segment (local system, whole chamber);Y, cm; entries",
01894                         nYbins,Ymin, Ymax));
01895             //
01896             Chan = 200;
01897             minChan = -0.2;
01898             maxChan = 0.2;
01899             sprintf(SpecName,"Phi_InefficientRecHits_inSegment_Ch%d_L%d",iChamber,iLayer);
01900             ChHist[ec][st][rg][iChamber-FirstCh].Phi_InefficientRecHits_inSegment.push_back
01901               (new TH1F(SpecName,"Missing RecHit/layer in a segment (local system, whole chamber);Phi, rad; entries",
01902                         Chan, minChan, maxChan));
01903             //
01904             sprintf(SpecName,"Phi_EfficientRecHits_inSegment_Ch%d_L%d",iChamber,iLayer);
01905             ChHist[ec][st][rg][iChamber-FirstCh].Phi_EfficientRecHits_inSegment.push_back
01906               (new TH1F(SpecName,"Efficient (extrapolated from the segment) in a segment (local system, whole chamber);Phi, rad; entries",
01907                         Chan, minChan, maxChan));
01908             
01909           }
01910           //
01911           sprintf(SpecName,"Sim_Rechits_Ch%d",iChamber);
01912           ChHist[ec][st][rg][iChamber-FirstCh].SimRechits = 
01913             new TH1F(SpecName,"Existing RecHit (Sim);layers (1-6);entries",nLayer_bins,Layer_min,Layer_max);
01914           //
01915           sprintf(SpecName,"Sim_Simhits_Ch%d",iChamber);
01916           ChHist[ec][st][rg][iChamber-FirstCh].SimSimhits = 
01917             new TH1F(SpecName,"Existing SimHit (Sim);layers (1-6);entries",nLayer_bins,Layer_min,Layer_max);
01918           //
01919           /*
01920           sprintf(SpecName,"Sim_Rechits_each_Ch%d",iChamber);
01921           ChHist[ec][st][rg][iChamber-FirstCh].SimRechits_each = 
01922             new TH1F(SpecName,"Existing RecHit (Sim), each;layers (1-6);entries",nLayer_bins,Layer_min,Layer_max);
01923           //
01924           sprintf(SpecName,"Sim_Simhits_each_Ch%d",iChamber);
01925           ChHist[ec][st][rg][iChamber-FirstCh].SimSimhits_each = 
01926             new TH1F(SpecName,"Existing SimHit (Sim), each;layers (1-6);entries",nLayer_bins,Layer_min,Layer_max);
01927           */
01928           theFile->cd();
01929         }
01930       }
01931     }
01932   }
01933 }
01934 
01935 // Destructor
01936 CSCEfficiency::~CSCEfficiency(){
01937   if (theService) delete theService; 
01938   // Write the histos to a file
01939   theFile->cd();
01940   //
01941   char SpecName[20];
01942   std::vector<float> bins, Efficiency, EffError;
01943   std::vector<float> eff(2);
01944 
01945   //---- loop over chambers
01946   std::map <std::string, bool> chamberTypes;
01947   chamberTypes["ME11"] = false;
01948   chamberTypes["ME12"] = false;
01949   chamberTypes["ME13"] = false;
01950   chamberTypes["ME21"] = false;
01951   chamberTypes["ME22"] = false;
01952   chamberTypes["ME31"] = false;
01953   chamberTypes["ME32"] = false;
01954   chamberTypes["ME41"] = false;
01955   
01956   map<std::string,bool>::iterator iter;  
01957   std::cout<<" Writing proper histogram structure (patience)..."<<std::endl;
01958   for(int ec = 0;ec<2;++ec){
01959     for(int st = 0;st<4;++st){
01960       sprintf(SpecName,"Stations__E%d_S%d",ec+1, st+1);
01961       theFile->cd(SpecName);
01962       StHist[ec][st].segmentChi2_ndf->Write();
01963       StHist[ec][st].hitsInSegment->Write();
01964       StHist[ec][st].AllSegments_eta->Write();
01965       StHist[ec][st].EfficientSegments_eta->Write();
01966       StHist[ec][st].ResidualSegments->Write();
01967       StHist[ec][st].EfficientSegments_XY->Write();
01968       StHist[ec][st].InefficientSegments_XY->Write();
01969       StHist[ec][st].EfficientALCT_momTheta->Write();
01970       StHist[ec][st].InefficientALCT_momTheta->Write();
01971       StHist[ec][st].EfficientCLCT_momPhi->Write();
01972       StHist[ec][st].InefficientCLCT_momPhi->Write();
01973       for(int rg = 0;rg<3;++rg){
01974         if(0!=st && rg>1){
01975           continue;
01976         }
01977         else if(1==rg && 3==st){
01978           continue;
01979         }
01980         for(int iChamber=FirstCh;iChamber<FirstCh+NumCh;iChamber++){
01981           if(0!=st && 0==rg && iChamber >18){
01982             continue;
01983           }
01984           sprintf(SpecName,"Chambers__E%d_S%d_R%d_Chamber_%d",ec+1, st+1, rg+1,iChamber);
01985           theFile->cd(SpecName);
01986           
01987           ChHist[ec][st][rg][iChamber-FirstCh].EfficientRechits_inSegment->Write();
01988           ChHist[ec][st][rg][iChamber-FirstCh].AllSingleHits->Write();
01989           ChHist[ec][st][rg][iChamber-FirstCh].digiAppearanceCount->Write();
01990           ChHist[ec][st][rg][iChamber-FirstCh].EfficientALCT_dydz->Write();
01991           ChHist[ec][st][rg][iChamber-FirstCh].InefficientALCT_dydz->Write();
01992           ChHist[ec][st][rg][iChamber-FirstCh].EfficientCLCT_dxdz->Write();
01993           ChHist[ec][st][rg][iChamber-FirstCh].InefficientCLCT_dxdz->Write();
01994           ChHist[ec][st][rg][iChamber-FirstCh].InefficientSingleHits->Write();
01995           ChHist[ec][st][rg][iChamber-FirstCh].EfficientRechits_good->Write();
01996           ChHist[ec][st][rg][iChamber-FirstCh].EfficientStrips->Write();
01997           ChHist[ec][st][rg][iChamber-FirstCh].StripWiresCorrelations->Write();
01998           ChHist[ec][st][rg][iChamber-FirstCh].NoWires_momTheta->Write();
01999           ChHist[ec][st][rg][iChamber-FirstCh].NoStrips_momPhi->Write();
02000           ChHist[ec][st][rg][iChamber-FirstCh].EfficientWireGroups->Write();
02001           for(unsigned int iLayer = 0; iLayer< 6; iLayer++){
02002             ChHist[ec][st][rg][iChamber-FirstCh].Y_InefficientRecHits_inSegment[iLayer]->Write();
02003             ChHist[ec][st][rg][iChamber-FirstCh].Y_EfficientRecHits_inSegment[iLayer]->Write();
02004             ChHist[ec][st][rg][iChamber-FirstCh].Phi_InefficientRecHits_inSegment[iLayer]->Write();
02005             ChHist[ec][st][rg][iChamber-FirstCh].Phi_EfficientRecHits_inSegment[iLayer]->Write();
02006           }
02007           ChHist[ec][st][rg][iChamber-FirstCh].SimRechits->Write();
02008           ChHist[ec][st][rg][iChamber-FirstCh].SimSimhits->Write();
02009           /*
02010           ChHist[ec][st][rg][iChamber-FirstCh].SimRechits_each->Write();
02011           ChHist[ec][st][rg][iChamber-FirstCh].SimSimhits_each->Write();
02012           */
02013           //
02014           theFile->cd(SpecName);
02015           theFile->cd();
02016         }
02017       }
02018     }
02019   }
02020   //
02021   sprintf(SpecName,"AllChambers");
02022   theFile->mkdir(SpecName);
02023   theFile->cd(SpecName);
02024   DataFlow->Write(); 
02025   TriggersFired->Write();
02026   ALCTPerEvent->Write();
02027   CLCTPerEvent->Write();
02028   recHitsPerEvent->Write();
02029   segmentsPerEvent->Write();
02030   //
02031   theFile->cd(SpecName);
02032   //---- Close the file
02033   theFile->Close();
02034 }
02035 
02036 // ------------ method called once each job just before starting event loop  ------------
02037 void
02038 CSCEfficiency::beginJob()
02039 {
02040 }
02041 
02042 // ------------ method called once each job just after ending the event loop  ------------
02043 void
02044 CSCEfficiency::endJob() {
02045 }
02046 
02047 DEFINE_FWK_MODULE(CSCEfficiency);
02048
Generated for CMSSW Reference Manual by doxygen 1.7.3