Reference Manual

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 //
00003 // Package:    SimHitShifter
00004 // Class:      SimHitShifter
00005 // 
00013 //
00014 // Original Author:  Camilo Andres Carrillo Montoya,40 2-B15,+41227671625,
00015 //         Created:  Mon Aug 30 18:35:05 CEST 2010
00016 // $Id: SimHitShifter.cc,v 1.1 2011/11/18 03:52:58 jiechen Exp $
00017 //
00018 //
00019 
00020 
00021 // system include files
00022 #include <memory>
00023 
00024 // user include files
00025 #include "FWCore/Framework/interface/Frameworkfwd.h"
00026 #include "FWCore/Framework/interface/EDProducer.h"
00027 
00028 #include "FWCore/Framework/interface/Event.h"
00029 #include "FWCore/Framework/interface/MakerMacros.h"
00030 
00031 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00032 
00033 #include "DataFormats/Common/interface/Handle.h"
00034 #include "FWCore/Framework/interface/ESHandle.h"
00035 #include "DataFormats/RPCDigi/interface/RPCDigi.h"
00036 #include "DataFormats/RPCDigi/interface/RPCDigiCollection.h"
00037 #include <DataFormats/RPCRecHit/interface/RPCRecHit.h>
00038 #include "DataFormats/RPCRecHit/interface/RPCRecHitCollection.h"
00039 #include <DataFormats/MuonDetId/interface/RPCDetId.h>
00040 #include <Geometry/RPCGeometry/interface/RPCGeometry.h>
00041 #include "Geometry/RPCGeometry/interface/RPCGeomServ.h"
00042 #include <DataFormats/GeometrySurface/interface/LocalError.h>
00043 #include <DataFormats/GeometryVector/interface/LocalPoint.h>
00044 #include "DataFormats/GeometrySurface/interface/Surface.h"
00045 #include "DataFormats/DetId/interface/DetId.h"
00046 #include <Geometry/Records/interface/MuonGeometryRecord.h>
00047 #include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
00048 #include "DataFormats/Candidate/interface/Candidate.h"
00049 #include "DataFormats/Candidate/interface/CandMatchMap.h"
00050 #include "DataFormats/Candidate/interface/CandidateFwd.h"
00051 #include "SimGeneral/HepPDTRecord/interface/ParticleDataTable.h"
00052 
00053 #include "MagneticField/Engine/interface/MagneticField.h"
00054 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
00055 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
00056 #include <DataFormats/GeometrySurface/interface/LocalError.h>
00057 #include <DataFormats/GeometryVector/interface/LocalPoint.h>
00058 #include "DataFormats/GeometrySurface/interface/Surface.h"
00059 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
00060 
00061 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
00062 #include "TrackingTools/GeomPropagators/interface/AnalyticalPropagator.h"
00063 #include "TrackPropagation/SteppingHelixPropagator/interface/SteppingHelixPropagator.h"
00064 
00065 #include "SimDataFormats/Track/interface/SimTrack.h"
00066 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
00067 #include "SimDataFormats/Vertex/interface/SimVertex.h"
00068 #include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
00069 #include "FastSimulation/Tracking/test/FastTrackAnalyzer.h"
00070 
00071 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
00072 #include "RecoMuon/TrackingTools/interface/MuonPatternRecoDumper.h"
00073 
00074 #include "DataFormats/MuonDetId/interface/MuonSubdetId.h"
00075 
00076 #include "DataFormats/MuonDetId/interface/CSCDetId.h"
00077 
00078 #include "Geometry/RPCGeometry/interface/RPCGeometry.h"
00079 #include <Geometry/RPCGeometry/interface/RPCRoll.h>
00080 #include <Geometry/Records/interface/MuonGeometryRecord.h>
00081 #include <Geometry/RPCGeometry/interface/RPCGeomServ.h>
00082 #include <Geometry/CommonDetUnit/interface/GeomDet.h>
00083 
00084 #include <cmath>
00085 
00086 //Root
00087 #include "TFile.h"
00088 #include "TF1.h"
00089 #include "TH1F.h"
00090 #include "TH1.h"
00091 #include "TH2F.h"
00092 #include "TROOT.h"
00093 #include "TMath.h"
00094 #include "TCanvas.h"
00095 
00096 //Track
00097 #include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h"
00098 #include "TrackingTools/PatternTools/interface/Trajectory.h"
00099 #include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h"
00100 #include "DataFormats/TrackCandidate/interface/TrackCandidate.h" 
00101 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
00102 
00103 #include<fstream>
00104 
00105 
00106 //
00107 // class declaration
00108 //
00109 
00110 class SimHitShifter : public edm::EDProducer {
00111    public:
00112       explicit SimHitShifter(const edm::ParameterSet&);
00113       ~SimHitShifter();
00114   //edm::ESHandle <RPCGeometry> rpcGeo;
00115       virtual void beginRun(const edm::Run&, const edm::EventSetup&);
00116       std::map<int,float> shiftinfo;
00117 
00118 
00119    private:
00120       std::string ShiftFileName;
00121       virtual void beginJob(const edm::Run&, const edm::EventSetup&) ;
00122       virtual void produce(edm::Event&, const edm::EventSetup&);
00123       virtual void endJob() ;
00124     
00125 };
00126 
00127 SimHitShifter::SimHitShifter(const edm::ParameterSet& iConfig)
00128 {
00129   std::cout<<"in the constructor"<<std::endl;
00130   
00131   ShiftFileName  = iConfig.getUntrackedParameter<std::string>("ShiftFileName","/afs/cern.ch/user/c/carrillo/simhits/CMSSW_3_5_8_patch2/src/simhitshifter/SimHitShifter/Merged_Muon_RawId_Shift.txt");
00132  
00133   //iSetup.get<MuonGeometryRecord>().get(rpcGeo);
00134 
00135   std::ifstream ifin(ShiftFileName.c_str());
00136 
00137   int rawId;
00138   float offset;
00139 
00140   std::cout<<"In the constructor, The name of the file is "<<ShiftFileName.c_str()<<std::endl;
00141 
00142   if(!ifin) std::cout<<"Problem reading the map rawId shift "<<ShiftFileName.c_str()<<std::endl;
00143   assert(ifin);
00144 
00145   while (ifin.good()){
00146     ifin >>rawId >>offset;
00147     shiftinfo[rawId]=offset;
00148     std::cout<<"rawId ="<<rawId<<" offset="<<offset<<std::endl;
00149   }
00150   
00151   produces<edm::PSimHitContainer>("MuonCSCHits");
00152   produces<edm::PSimHitContainer>("MuonDTHits");
00153   produces<edm::PSimHitContainer>("MuonRPCHits");
00154 }
00155 
00156 
00157 SimHitShifter::~SimHitShifter()
00158 {
00159 }
00160 
00161 void SimHitShifter::produce(edm::Event& iEvent, const edm::EventSetup& iSetup){
00162    using namespace edm;
00163 
00164    //std::cout << " Getting the SimHits " <<std::endl;
00165    std::vector<edm::Handle<edm::PSimHitContainer> > theSimHitContainers;
00166    iEvent.getManyByType(theSimHitContainers);
00167    //std::cout << " The Number of sim Hits is  " << theSimHitContainers.size() <<std::endl;
00168 
00169    std::auto_ptr<edm::PSimHitContainer> pcsc(new edm::PSimHitContainer);
00170    std::auto_ptr<edm::PSimHitContainer> pdt(new edm::PSimHitContainer);
00171    std::auto_ptr<edm::PSimHitContainer> prpc(new edm::PSimHitContainer);
00172 
00173    std::vector<PSimHit> theSimHits;
00174 
00175    using std::oct;
00176    using std::dec;
00177    
00178    for (int i = 0; i < int(theSimHitContainers.size()); i++){
00179      theSimHits.insert(theSimHits.end(),theSimHitContainers.at(i)->begin(),theSimHitContainers.at(i)->end());
00180    } 
00181 
00182    for (std::vector<PSimHit>::const_iterator iHit = theSimHits.begin(); iHit != theSimHits.end(); iHit++){
00183      DetId theDetUnitId((*iHit).detUnitId());
00184      DetId simdetid= DetId((*iHit).detUnitId());
00185 
00186      if(simdetid.det()!=DetId::Muon) continue;
00187 
00188      float newtof = 0;
00189     
00190      if(simdetid.det()==DetId::Muon &&  simdetid.subdetId()== MuonSubdetId::RPC){//Only RPCs
00191        //std::cout<<"\t\t We have an RPC Sim Hit! in t="<<(*iHit).timeOfFlight()<<" DetId="<<(*iHit).detUnitId()<<std::endl;
00192        if(shiftinfo.find(simdetid.rawId())==shiftinfo.end()){
00193          std::cout<<"RPC Warning the RawId = "<<simdetid.det()<<" | "<<simdetid.rawId()<<"is not in the map"<<std::endl;
00194          newtof = (*iHit).timeOfFlight();
00195        }else{
00196          newtof = (*iHit).timeOfFlight()+shiftinfo[simdetid.rawId()];
00197        }
00198        
00199        PSimHit hit((*iHit).entryPoint(),(*iHit).exitPoint(),(*iHit).pabs(),
00200                    newtof,
00201                    (*iHit).energyLoss(),(*iHit).particleType(),simdetid,(*iHit). trackId(),(*iHit).thetaAtEntry(),(*iHit).phiAtEntry(),(*iHit).processType());
00202        prpc->push_back(hit);
00203      }
00204      else if(simdetid.det()==DetId::Muon &&  simdetid.subdetId()== MuonSubdetId::DT){//Only DTs
00205        int RawId = simdetid.rawId(); 
00206        std::cout<<"We found a DT simhit the RawId in Dec is";
00207        std::cout<<dec<<RawId<<std::endl;
00208        std::cout<<"and in oct"<<std::endl;
00209        std::cout<<oct<<RawId<< std::endl;
00210        std::cout<<"once masked in oct "<<std::endl;
00211        int compressedRawId = RawId/8/8/8/8/8;
00212        std::cout<<compressedRawId<<std::endl;
00213        std::cout<<"extendedRawId"<<std::endl;
00214        int extendedRawId = compressedRawId*8*8*8*8*8;
00215        std::cout<<extendedRawId<<std::endl;
00216        std::cout<<"converted again in decimal"<<std::endl;
00217        std::cout<<dec<<extendedRawId<<std::endl;
00218        
00219        if(shiftinfo.find(extendedRawId)==shiftinfo.end()){
00220          //std::cout<<"DT Warning the RawId = "<<extendedRawId<<"is not in the map"<<std::endl;
00221          newtof = (*iHit).timeOfFlight();
00222        }else{
00223          newtof = (*iHit).timeOfFlight()+shiftinfo[extendedRawId];
00224          std::cout<<"RawId = "<<extendedRawId<<"is in the map "<<(*iHit).timeOfFlight()<<" "<<newtof<<std::endl;
00225        }
00226        
00227        std::cout<<"\t\t We have an DT Sim Hit! in t="<<(*iHit).timeOfFlight()<<" DetId="<<(*iHit).detUnitId()<<std::endl;
00228        PSimHit hit((*iHit).entryPoint(),(*iHit).exitPoint(),(*iHit).pabs(),
00229                    newtof,
00230                    (*iHit).energyLoss(),(*iHit).particleType(),simdetid,(*iHit). trackId(),(*iHit).thetaAtEntry(),(*iHit).phiAtEntry(),(*iHit).processType());
00231        pdt->push_back(hit);
00232      }
00233      else if(simdetid.det()==DetId::Muon &&  simdetid.subdetId()== MuonSubdetId::CSC){//Only CSCs
00234        //std::cout<<"\t\t We have an CSC Sim Hit! in t="<<(*iHit).timeOfFlight()<<" DetId="<<(*iHit).detUnitId()<<std::endl;
00235        
00236        CSCDetId TheCSCDetId = CSCDetId(simdetid);
00237        CSCDetId TheChamberDetId = TheCSCDetId.chamberId();
00238        
00239        if(shiftinfo.find(TheChamberDetId.rawId())==shiftinfo.end()){
00240          std::cout<<"The RawId is not in the map,perhaps it is on the CSCs station 1 ring 4"<<std::endl;
00241          if(TheChamberDetId.station()==1 && TheChamberDetId.ring()==4){
00242            CSCDetId TheChamberDetIdNoring4= CSCDetId(TheChamberDetId.endcap(),TheChamberDetId.station(),1 //1 instead of 4
00243                                                      ,TheChamberDetId.chamber(),TheChamberDetId.layer());
00244            
00245            if(shiftinfo.find(TheChamberDetIdNoring4.rawId())==shiftinfo.end()){
00246              std::cout<<"CSC Warning the RawId = "<<TheChamberDetIdNoring4<<" "<<TheChamberDetIdNoring4.rawId()<<"is not in the map"<<std::endl;
00247              newtof = (*iHit).timeOfFlight();
00248            }else{
00249              newtof = (*iHit).timeOfFlight()+shiftinfo[TheChamberDetIdNoring4.rawId()];
00250            }
00251          }
00252        }else{
00253          newtof = (*iHit).timeOfFlight()+shiftinfo[TheChamberDetId.rawId()];
00254        }
00255        
00256        PSimHit hit((*iHit).entryPoint(),(*iHit).exitPoint(),(*iHit).pabs(),
00257                    newtof,
00258                    (*iHit).energyLoss(),(*iHit).particleType(),simdetid,(*iHit). trackId(),(*iHit).thetaAtEntry(),(*iHit).phiAtEntry(),(*iHit).processType());
00259        
00260        std::cout<<"CSC check newtof"<<newtof<<" "<<(*iHit).timeOfFlight()<<std::endl;
00261        if(newtof==(*iHit).timeOfFlight())std::cout<<"Warning!!!"<<std::endl;
00262        pcsc->push_back(hit);
00263      }     
00264    }
00265 
00266    std::cout<<"Putting collections in the event"<<std::endl;
00267 
00268    iEvent.put(pcsc,"MuonCSCHits");
00269    iEvent.put(pdt,"MuonDTHits");
00270    iEvent.put(prpc,"MuonRPCHits");
00271    
00272 }
00273 
00274 void 
00275 SimHitShifter::beginRun(const edm::Run& run, const edm::EventSetup& iSetup)
00276 {
00277 
00278 }
00279 
00280 // ------------ method called once each job just before starting event loop  ------------
00281 void 
00282 SimHitShifter::beginJob(const edm::Run& run, const edm::EventSetup& iSetup)
00283 {
00284 
00285 }
00286 
00287 // ------------ method called once each job just after ending the event loop  ------------
00288 void 
00289 SimHitShifter::endJob() {
00290 }
00291 
00292 //define this as a plug-in
00293 DEFINE_FWK_MODULE(SimHitShifter);
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_5_3_14/src/SUSYBSMAnalysis/HSCP/src/SimHitShifter.cc
