From 8387c4fea7c83bea333a04cf8b394efa2412b92f Mon Sep 17 00:00:00 2001 From: Alessio Bonato <sha1-45eb77c4b8d9753c44d0fc05483b6faa1620d239@cern.ch> Date: Fri, 2 Oct 2009 12:56:09 +0000 Subject: [PATCH] --- yaml --- r: 74603 b: "refs/heads/CMSSW_7_1_X" c: 8c189c34b093970176470cc9835afc9b3fce0931 h: "refs/heads/CMSSW_7_1_X" i: 74601: 0eb0edf105cd05610fe6afd31825ead860530be4 74599: 574c7788a7558fe0e41cd6f8dd4fff05a7548791 v: v3 --- [refs] | 2 +- trunk/Alignment/TrackerAlignment/BuildFile | 8 +- .../plugins/AlignmentPrescaler.cc | 242 ++++++++++++++++ .../plugins/AlignmentPrescaler.h | 80 ++++++ .../TrackerAlignment/plugins/BuildFile | 28 +- .../TrackerAlignment/plugins/OverlapTagger.cc | 204 ++++++++++++++ .../TrackerAlignment/plugins/OverlapTagger.h | 51 ++++ .../TrackerAlignment/plugins/TreeMerger.cc | 260 ++++++++++++++++++ .../TrackerAlignment/plugins/TreeMerger.h | 59 ++++ .../python/AlignmentPrescaler_cff.py | 8 + .../python/OverlapTagger_cff.py | 9 + .../TrackerAlignment/python/TreeMerger_cff.py | 16 ++ 12 files changed, 961 insertions(+), 6 deletions(-) create mode 100644 trunk/Alignment/TrackerAlignment/plugins/AlignmentPrescaler.cc create mode 100644 trunk/Alignment/TrackerAlignment/plugins/AlignmentPrescaler.h create mode 100644 trunk/Alignment/TrackerAlignment/plugins/OverlapTagger.cc create mode 100644 trunk/Alignment/TrackerAlignment/plugins/OverlapTagger.h create mode 100644 trunk/Alignment/TrackerAlignment/plugins/TreeMerger.cc create mode 100644 trunk/Alignment/TrackerAlignment/plugins/TreeMerger.h create mode 100644 trunk/Alignment/TrackerAlignment/python/AlignmentPrescaler_cff.py create mode 100644 trunk/Alignment/TrackerAlignment/python/OverlapTagger_cff.py create mode 100644 trunk/Alignment/TrackerAlignment/python/TreeMerger_cff.py diff --git a/[refs] b/[refs] index e48dc912a33d8..542247c9abb71 100644 --- a/[refs] +++ b/[refs] @@ -1,3 +1,3 @@ --- refs/heads/gh-pages: 09c786f70121f131b3715aaf3464996502bbeb7e -"refs/heads/CMSSW_7_1_X": e5c23c20f8f591c417c106e295d101d907840b26 +"refs/heads/CMSSW_7_1_X": 8c189c34b093970176470cc9835afc9b3fce0931 diff --git a/trunk/Alignment/TrackerAlignment/BuildFile b/trunk/Alignment/TrackerAlignment/BuildFile index eb0fc32515bd0..d56c05970bcff 100644 --- a/trunk/Alignment/TrackerAlignment/BuildFile +++ b/trunk/Alignment/TrackerAlignment/BuildFile @@ -13,8 +13,11 @@ <use name=Geometry/CommonDetUnit> <use name=Geometry/Records> <use name=Geometry/TrackerGeometryBuilder> +<use name=DataFormats/Alignment> +<use name=AnalysisDataFormats/SiStripClusterInfo> +<use name=RecoTracker/TransientTrackingRecHit> <export> - <use name=Alignment/CommonAlignment> +# <use name=Alignment/CommonAlignment> <use name=CondFormats/Alignment> <use name=DataFormats/SiPixelDetId> <use name=DataFormats/SiStripDetId> @@ -22,6 +25,3 @@ <lib name=AlignmentTrackerAlignment> </export> - - - diff --git a/trunk/Alignment/TrackerAlignment/plugins/AlignmentPrescaler.cc b/trunk/Alignment/TrackerAlignment/plugins/AlignmentPrescaler.cc new file mode 100644 index 0000000000000..cdb58f5e0b5f3 --- /dev/null +++ b/trunk/Alignment/TrackerAlignment/plugins/AlignmentPrescaler.cc @@ -0,0 +1,242 @@ +#include "Alignment/TrackerAlignment/plugins/AlignmentPrescaler.h" + + +AlignmentPrescaler::AlignmentPrescaler(const edm::ParameterSet &iConfig): + src_(iConfig.getParameter<edm::InputTag>("src")), + AM_(iConfig.getParameter<edm::InputTag>("assomap")), + prescfilename_(iConfig.getParameter<std::string>("PrescFileName")), + presctreename_(iConfig.getParameter<std::string>("PrescTreeName")) +{ + // issue the produce<> + produces<AliClusterValueMap>(); + produces<AliTrackTakenClusterValueMap>(); + +} + +AlignmentPrescaler::~AlignmentPrescaler(){ + // +} + +void AlignmentPrescaler::beginJob( const edm::EventSetup & ){ + // + fpresc_=new TFile(prescfilename_.c_str(),"READ"); + tpresc_=(TTree*)fpresc_->Get(presctreename_.c_str()); + tpresc_->BuildIndex("DetId"); + tpresc_->SetBranchStatus("*",0); + tpresc_->SetBranchStatus("DetId",1); + tpresc_->SetBranchStatus("PrescaleFactor",1); + tpresc_->SetBranchStatus("PrescaleFactorOverlap",1); + + detid_=0; + preschit_=99.0; + prescoverlap_=88.0; + + tpresc_->SetBranchAddress("DetId",&detid_); + tpresc_->SetBranchAddress("PrescaleFactor",&preschit_); + tpresc_->SetBranchAddress("PrescaleFactorOverlap",&prescoverlap_); + + myrand=new TRandom3(); + // myrand->SetSeed(); + + + /* + fh=new TFile("debugfile1.root","RECREATE"); + hrr=new TH1F("randdistr","distribution of random numbers",100,0.0,1.0); + hnhitssubdet=new TH1F("hitvssubdet","Distribution of hits Vs SubDet",7,0.0,7.0); + hnhitsTIBL3=new TH1I("hitvsphiTIBL3","Distribution of hits Vs DetId (TIB L3 only)",5848,369153044,369158892); + + totnhitspxl_=0; + */ + + +} + +void AlignmentPrescaler::endJob( ){ + /* + cout<<"\n\n%%%%%%%% At the end of AlignmentPrescale the number of PIXEL HITS TAKEN was "<< totnhitspxl_<<endl<<endl; + + fh->cd(); + hrr->Write(); + hnhitssubdet->Write(); + hnhitsTIBL3->Write(); + + + delete hrr; + fh->Close(); + delete fh; + */ + + // + delete tpresc_; + fpresc_->Close(); + delete fpresc_; + delete myrand; +} + +void AlignmentPrescaler::produce(edm::Event &iEvent, const edm::EventSetup &iSetup){ + // std::cout<<"\n\n#################\n### Starting the AlignmentPrescaler::produce ; Event: "<<iEvent.id().run() <<", "<<iEvent.id().event()<<std::endl; + edm::Handle<reco::TrackCollection> Tracks; + iEvent.getByLabel(src_, Tracks); + + //take HitAssomap + edm::Handle<AliClusterValueMap> hMap; + iEvent.getByLabel(AM_, hMap); + AliClusterValueMap InValMap=*hMap; + + //prepare the output of the ValueMap flagging tracks + std::vector<int> trackflags(Tracks->size(),0); + + + //int npxlhits=0; + + //loop on tracks + for(std::vector<reco::Track>::const_iterator ittrk = Tracks->begin(), edtrk = Tracks->end(); ittrk != edtrk; ++ittrk){ + //loop on tracking rechits + // std::cout << "Loop on hits of track #" << (ittrk - Tracks->begin()) << std::endl; + int nhit=0; + int ntakenhits=0; + bool firstTakenHit=false; + + for (trackingRecHit_iterator ith = ittrk->recHitsBegin(), edh = ittrk->recHitsEnd(); ith != edh; ++ith) { + const TrackingRecHit *hit = ith->get(); // ith is an iterator on edm::Ref to rechit + if(! hit->isValid()){ + nhit++; + continue; + } + uint32_t tmpdetid = hit->geographicalId().rawId(); + tpresc_->GetEntryWithIndex(tmpdetid); + + + //------------- + //decide whether to take this hit or not + bool takeit=false; + int subdetId=hit->geographicalId().subdetId(); + + + //check first if the cluster is also in the overlap asso map + bool isOverlapHit=false; + // bool first=true; + //ugly... + const SiStripRecHit2D* striphit = dynamic_cast<const SiStripRecHit2D*>(hit); + const SiPixelRecHit* pixelhit= dynamic_cast<const SiPixelRecHit*>(hit); + AlignmentClusterFlag tmpflag(hit->geographicalId()); + if(subdetId>2){// SST case + if(striphit!=0){ + SiStripRecHit2D::ClusterRef stripclust(striphit->cluster()); + tmpflag=InValMap[stripclust]; + tmpflag.SetDetId(hit->geographicalId()); + if(tmpflag.isOverlap())isOverlapHit=true; + // cout<<"~*~*~* Prescale for module "<<tmpflag.detId().rawId()<<"("<<InValMap[stripclust].detId().rawId() <<") is "<<preschit_<<flush; + //if(tmpflag.isOverlap())cout<<" (it is Overlap)"<<flush; + // else cout<<endl; + + }//end if striphit!=0 + }//end if is a strip hit + else{ + if(pixelhit!=0){ + //npxlhits++; + SiPixelClusterRefNew pixclust(pixelhit->cluster()); + tmpflag=InValMap[pixclust]; + tmpflag.SetDetId(hit->geographicalId()); + if(tmpflag.isOverlap())isOverlapHit=true; + } + }//end else is a pixel hit + // tmpflag.SetDetId(hit->geographicalId()); + + if( isOverlapHit ){ + //cout<<" DetId="<<tmpdetid<<" is Overlap! "<<flush; + takeit=(float(myrand->Rndm())<=prescoverlap_); + } + if( !takeit ){ + float rr=float(myrand->Rndm()); + takeit=(rr<=preschit_); + } + if(takeit){//HIT TAKEN ! + //cout<<" DetId="<<tmpdetid<<" taken!"<<flush; + tmpflag.SetTakenFlag(); + + if(subdetId>2){ + SiStripRecHit2D::ClusterRef stripclust(striphit->cluster()); + InValMap[stripclust]=tmpflag;//.SetTakenFlag(); + + } + else{ + SiPixelClusterRefNew pixclust(pixelhit->cluster()); + InValMap[pixclust]=tmpflag;//.SetTakenFlag(); + } + + if(!firstTakenHit){ + firstTakenHit=true; + //std::cout<<"Index of the track iterator is "<< ittrk-Tracks->begin() <<endl; + + } + ntakenhits++; + }//end if take this hit + //cout<<endl; + + nhit++; + //cout<<endl; + }//end loop on RecHits + trackflags[ittrk-Tracks->begin()]=ntakenhits; + // cout<<"Entrioes in debug histo: "<<hrr->GetEntries()<<endl; + }//end loop on tracks + + + + // totnhitspxl_+=ntakenhits; + //cout<<"AlignmentPrescaler::produce says that in this event "<<ntakenhits<<" pixel clusters were taken (out of "<<npxlhits<<" total pixel hits."<<endl; + + + + //save the asso map, tracks... + // prepare output + std::auto_ptr<AliClusterValueMap> OutVM( new AliClusterValueMap); + *OutVM=InValMap; + + iEvent.put(OutVM); + + + std::auto_ptr<AliTrackTakenClusterValueMap> trkVM( new AliTrackTakenClusterValueMap); + AliTrackTakenClusterValueMap::Filler trkmapfiller(*trkVM); + trkmapfiller.insert(Tracks,trackflags.begin(),trackflags.end() ); + trkmapfiller.fill(); + iEvent.put(trkVM); + + +}//end produce + + +int AlignmentPrescaler::layerFromId (const DetId& id) const +{ + if ( uint32_t(id.subdetId())==PixelSubdetector::PixelBarrel ) { + PXBDetId tobId(id); + return tobId.layer(); + } + else if ( uint32_t(id.subdetId())==PixelSubdetector::PixelEndcap ) { + PXFDetId tobId(id); + return tobId.disk() + (3*(tobId.side()-1)); + } + else if ( id.subdetId()==StripSubdetector::TIB ) { + TIBDetId tibId(id); + return tibId.layer(); + } + else if ( id.subdetId()==StripSubdetector::TOB ) { + TOBDetId tobId(id); + return tobId.layer(); + } + else if ( id.subdetId()==StripSubdetector::TEC ) { + TECDetId tobId(id); + return tobId.wheel() + (9*(tobId.side()-1)); + } + else if ( id.subdetId()==StripSubdetector::TID ) { + TIDDetId tobId(id); + return tobId.wheel() + (3*(tobId.side()-1)); + } + return -1; + +}//end layerfromId + +// ========= MODULE DEF ============== +#include "FWCore/PluginManager/interface/ModuleDef.h" +#include "FWCore/Framework/interface/MakerMacros.h" +DEFINE_FWK_MODULE(AlignmentPrescaler); diff --git a/trunk/Alignment/TrackerAlignment/plugins/AlignmentPrescaler.h b/trunk/Alignment/TrackerAlignment/plugins/AlignmentPrescaler.h new file mode 100644 index 0000000000000..3e50bd4ecffcb --- /dev/null +++ b/trunk/Alignment/TrackerAlignment/plugins/AlignmentPrescaler.h @@ -0,0 +1,80 @@ +#ifndef TrackerAlignment_AlignmentPrescaler_H +#define TrackerAlignment_AlignmentPrescaler_H + +#include <Riostream.h> +#include <string> +#include "TFile.h" +#include "TTree.h" +#include "TRandom3.h" +#include "TH1F.h" + + +#include "FWCore/Framework/interface/EDProducer.h" +#include "FWCore/Framework/interface/EventPrincipal.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/EventSetup.h" +#include "FWCore/Framework/interface/ESHandle.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/ParameterSet/interface/InputTag.h" + +#include "DataFormats/Common/interface/View.h" +#include "DataFormats/DetId/interface/DetId.h" +#include "DataFormats/SiPixelDetId/interface/PXBDetId.h" +#include "DataFormats/SiPixelDetId/interface/PXFDetId.h" +#include "DataFormats/SiStripDetId/interface/TIBDetId.h" +#include "DataFormats/SiStripDetId/interface/TIDDetId.h" +#include "DataFormats/SiStripDetId/interface/TOBDetId.h" +#include "DataFormats/SiStripDetId/interface/TECDetId.h" +#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h" +#include "Alignment/TrackerAlignment/interface/AlignableTracker.h" +#include "Alignment/CommonAlignment/interface/Alignable.h" +#include "Alignment/CommonAlignment/interface/Utilities.h" +#include "Utilities/General/interface/ClassName.h" + +#include "DataFormats/TrackReco/interface/Track.h" +#include "AnalysisDataFormats/TrackInfo/interface/TrackInfo.h" + +#include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2D.h" +#include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHit.h" +#include "RecoTracker/TransientTrackingRecHit/interface/TSiStripRecHit2DLocalPos.h" +#include "RecoTracker/TransientTrackingRecHit/interface/TSiPixelRecHit.h" +#include "DataFormats/Alignment/interface/AlignmentClusterFlag.h" +#include "DataFormats/Alignment/interface/AliClusterValueMap.h" + + +class AlignmentPrescaler : public edm::EDProducer{ + + public: + AlignmentPrescaler(const edm::ParameterSet &iConfig); + ~AlignmentPrescaler(); + void beginJob( const edm::EventSetup & ); + void endJob(); + virtual void produce(edm::Event &iEvent, const edm::EventSetup &iSetup) ; + + private: + edm::InputTag src_;//tracks in input + edm::InputTag AM_;//Hit-quality association map + + std::string prescfilename_;//name of the file containing the TTree with the prescaling factors + std::string presctreename_;//name of the TTree with the prescaling factors + + TFile *fpresc_; + TTree *tpresc_; + TRandom3 *myrand; + /* + //temp, just for debug + TFile *fh; + TH1F *hrr; + TH1F *hnhitssubdet; + TH1I *hnhitsTIBL3; + */ + + + int layerFromId (const DetId& id) const; + + unsigned int detid_; + float preschit_, prescoverlap_; + int totnhitspxl_; +}; +#endif diff --git a/trunk/Alignment/TrackerAlignment/plugins/BuildFile b/trunk/Alignment/TrackerAlignment/plugins/BuildFile index 45fff64c3d69e..f2a3fdfbb3076 100644 --- a/trunk/Alignment/TrackerAlignment/plugins/BuildFile +++ b/trunk/Alignment/TrackerAlignment/plugins/BuildFile @@ -5,4 +5,30 @@ <use name=Geometry/TrackingGeometryAligner> <use name=CondCore/DBOutputService> -<library name=AlignmentTrackerAlignmentPlugin file="*.cc"><flags EDM_PLUGIN=1></library> +<library name=AlignmentTrackerAlignmentPlugin file="*Misalign*.cc"> +<flags EDM_PLUGIN=1> +</library> + + + +#<use name=FWCore/MessageLogger> +#<use name=FWCore/Framework> +#<use name=FWCore/ParameterSet> + +<library name=AlignmentPrescalerPlugin file=AlignmentPrescaler.cc> + <use name=DataFormats/TrackerRecHit2D> + <use name=DataFormats/Alignment> + <flags EDM_PLUGIN=1> +</library> + +<library name="TreeMergerPlugin" file="TreeMerger.cc"> +# <use name=Alignment/TrackerAlignment> + <flags EDM_PLUGIN=1> +</library> + +<library name=OverlapTaggerPlugin file=OverlapTagger.cc> + <use name=DataFormats/TrackReco> + <use name=DataFormats/TrackerRecHit2D> + <use name=DataFormats/Alignment> + <flags EDM_PLUGIN=1> +</library> diff --git a/trunk/Alignment/TrackerAlignment/plugins/OverlapTagger.cc b/trunk/Alignment/TrackerAlignment/plugins/OverlapTagger.cc new file mode 100644 index 0000000000000..16c3616e5ce89 --- /dev/null +++ b/trunk/Alignment/TrackerAlignment/plugins/OverlapTagger.cc @@ -0,0 +1,204 @@ +#include "Alignment/TrackerAlignment/plugins/OverlapTagger.h" + +using namespace edm; +using namespace reco; + +OverlapTagger::OverlapTagger(const edm::ParameterSet& iConfig): + src_( iConfig.getParameter<edm::InputTag>("src") ), + srcClust_( iConfig.getParameter<edm::InputTag>("Clustersrc") ), + rejectBadMods_( iConfig.getParameter<bool>("rejectBadMods")), + BadModsList_( iConfig.getParameter<std::vector<uint32_t> >("BadMods")) +{ + + produces<AliClusterValueMap>(); //produces the ValueMap (VM) to be stored in the Event at the end +} + +OverlapTagger::~OverlapTagger(){} + + +void OverlapTagger::produce(edm::Event &iEvent, const edm::EventSetup &iSetup){ + edm::Handle<TrajTrackAssociationCollection> assoMap; + iEvent.getByLabel(src_, assoMap); + // cout<<"\n\n############################\n### Starting a new OverlapTagger - Ev "<<iEvent.id().run()<<", "<<iEvent.id().event()<<endl; + + AlignmentClusterFlag iniflag; + edm::Handle<edmNew::DetSetVector<SiPixelCluster> > pixelclusters; + iEvent.getByLabel(srcClust_, pixelclusters);//same label as tracks + std::vector<AlignmentClusterFlag> pixelvalues(pixelclusters->dataSize(), iniflag);//vector where to store value to be fileld in the VM + + edm::Handle<edmNew::DetSetVector<SiStripCluster> > stripclusters; + iEvent.getByLabel(srcClust_, stripclusters);//same label as tracks + std::vector<AlignmentClusterFlag> stripvalues(stripclusters->dataSize(), iniflag);//vector where to store value to be fileld in the VM + + + + //start doing the thing! + + //loop over trajectories + for (TrajTrackAssociationCollection::const_iterator itass = assoMap->begin(); itass != assoMap->end(); ++itass){ + + int nOverlaps=0; + const edm::Ref<std::vector<Trajectory> >traj = itass->key;//trajectory in the collection + const Trajectory * myTrajectory= &(*traj); + std::vector<TrajectoryMeasurement> tmColl =myTrajectory->measurements(); + + const reco::TrackRef tkref = itass->val;//associated track track in the collection + // const Track * trk = &(*tkref); + int hitcnt=-1; + + //loop over traj meas + const TrajectoryMeasurement* previousTM(0); + DetId previousId(0); + int previousLayer(-1); + + for(std::vector<TrajectoryMeasurement>::const_iterator itTrajMeas = tmColl.begin(); itTrajMeas!=tmColl.end(); ++itTrajMeas){ + hitcnt++; + + if ( previousTM!=0 ) { + // std::cout<<"Checking TrajMeas ("<<hitcnt+1<<"):"<<std::endl; + if(!previousTM->recHit()->isValid()){ + //std::cout<<"Previous RecHit invalid !"<<std::endl; + continue;} + //else std::cout<<"\nDetId: "<<std::flush<<previousTM->recHit()->geographicalId().rawId()<<"\t Local x of hit: "<<previousTM->recHit()->localPosition().x()<<std::endl; + } + else{ + //std::cout<<"This is the first Traj Meas of the Trajectory! The Trajectory contains "<< tmColl.size()<<" TrajMeas"<<std::endl; + } + + + TransientTrackingRecHit::ConstRecHitPointer hitpointer = itTrajMeas->recHit(); + const TrackingRecHit *hit=&(* hitpointer); + if(!hit->isValid())continue; + + //std::cout << " hit number " << (ith - itt->recHitsBegin()) << std::endl; + DetId detid = hit->geographicalId(); + int layer(layerFromId(detid));//layer 1-4=TIB, layer 5-10=TOB + int subDet = detid.subdetId(); + + + if ( ( previousTM!=0 )&& (layer!=-1 )) { + for (std::vector<TrajectoryMeasurement>::const_iterator itmCompare =itTrajMeas-1;itmCompare >= tmColl.begin() && itmCompare > itTrajMeas - 4;--itmCompare){ + DetId compareId = itmCompare->recHit()->geographicalId(); + if ( subDet != compareId.subdetId() || layer != layerFromId(compareId)) break; + if (!itmCompare->recHit()->isValid()) continue; + if ( (subDet<=2) || (subDet > 2 && SiStripDetId(detid).stereo()==SiStripDetId(compareId).stereo())){//if either pixel or strip stereo module + + // + //WOW, we have an overlap!!!!!! + // + AlignmentClusterFlag hitflag(hit->geographicalId()); + hitflag.SetOverlapFlag(); + // cout<<"Overlap found in SubDet "<<subDet<<"!!!"<<flush; + + if (subDet>2){ + //cout<<" TypeId of the RecHit: "<<className(*hit)<<endl; + const TSiStripRecHit2DLocalPos* transstriphit = dynamic_cast<const TSiStripRecHit2DLocalPos*>(hit); + + if(transstriphit!=0){ + const SiStripRecHit2D* striphit=transstriphit->specificHit(); + // cout<<"Pointer to SiStripRecHit2D= "<<striphit<<endl; + SiStripRecHit2D::ClusterRef stripclust(striphit->cluster()); + + if(stripclust.id()==stripclusters.id()){//ensure that the pixclust is really present in the original cluster collection!!! + stripvalues[stripclust.key()]=hitflag; + + //cout<<">>> Storing in the ValueMap a StripClusterRef with Cluster.Key: "<<stripclust.key()<<" ("<<striphit->cluster().key() <<"), Cluster.Id: "<<stripclust.id()<<" (DetId is "<<hit->geographicalId().rawId()<<")"<<endl; + } + else{ + cout<<"ERROR in <OverlapTagger::produce>: ProdId of Strip clusters mismatched: "<<stripclust.id()<<" vs "<<stripclusters.id() <<endl; + } + + // cout<<"Cluster baricentre: "<<stripclust->barycenter()<<endl; + } + else{ + cout<<"ERROR in <OverlapTagger::produce>: Dynamic cast of Strip RecHit failed! TypeId of the RecHit: "<<className(*hit)<<endl; + } + } + else {//pixel hit + const TSiPixelRecHit* transpixelhit = dynamic_cast<const TSiPixelRecHit*>(hit); + if(transpixelhit!=0){ + const SiPixelRecHit* pixelhit=transpixelhit->specificHit(); + SiPixelClusterRefNew pixclust(pixelhit->cluster()); + + if(pixclust.id()==pixelclusters.id()){ + pixelvalues[pixclust.key()]=hitflag; + //cout<<">>> Storing in the ValueMap a PixelClusterRef with ProdID: "<<pixclust.id()<<" (DetId is "<<hit->geographicalId().rawId()<<")" <<endl;//" and a Val with ID: "<<flag.id()<<endl; + } + else{ + cout<<"ERROR in <OverlapTagger::produce>: ProdId of Pixel clusters mismatched: "<<pixclust.id()<<" vs "<<pixelclusters.id() <<endl; + } + } + else{ + cout<<"ERROR in <OverlapTagger::produce>: Dynamic cast of Pixel RecHit failed! TypeId of the RecHit: "<<className(*hit)<<endl; + } + }//end 'else' it is a pixel hit + + nOverlaps++; + break; + } + }//end second loop on TM + }//end if a previous TM exists + + previousTM = &(* itTrajMeas); + previousId = detid; + previousLayer = layer; + }//end loop over traj meas + //std::cout<<"Found "<<nOverlaps<<" overlaps in this trajectory"<<std::endl; + + }//end loop over trajectories + + + + // prepare output + std::auto_ptr<AliClusterValueMap> hitvalmap( new AliClusterValueMap); + AliClusterValueMap::Filler mapfiller(*hitvalmap); + + edm::TestHandle<std::vector<AlignmentClusterFlag> > fakePixelHandle( &pixelvalues,pixelclusters.id()); + mapfiller.insert(fakePixelHandle, pixelvalues.begin(), pixelvalues.end()); + + edm::TestHandle<std::vector<AlignmentClusterFlag> > fakeStripHandle( &stripvalues,stripclusters.id()); + mapfiller.insert(fakeStripHandle, stripvalues.begin(), stripvalues.end()); + mapfiller.fill(); + + + + + + // iEvent.put(stripmap); + iEvent.put(hitvalmap); +}//end OverlapTagger::produce +int OverlapTagger::layerFromId (const DetId& id) const +{ + if ( uint32_t(id.subdetId())==PixelSubdetector::PixelBarrel ) { + PXBDetId tobId(id); + return tobId.layer(); + } + else if ( uint32_t(id.subdetId())==PixelSubdetector::PixelEndcap ) { + PXFDetId tobId(id); + return tobId.disk() + (3*(tobId.side()-1)); + } + else if ( id.subdetId()==StripSubdetector::TIB ) { + TIBDetId tibId(id); + return tibId.layer(); + } + else if ( id.subdetId()==StripSubdetector::TOB ) { + TOBDetId tobId(id); + return tobId.layer(); + } + else if ( id.subdetId()==StripSubdetector::TEC ) { + TECDetId tobId(id); + return tobId.wheel() + (9*(tobId.side()-1)); + } + else if ( id.subdetId()==StripSubdetector::TID ) { + TIDDetId tobId(id); + return tobId.wheel() + (3*(tobId.side()-1)); + } + return -1; + +}//end layerfromId + +// ========= MODULE DEF ============== +#include "FWCore/PluginManager/interface/ModuleDef.h" +#include "FWCore/Framework/interface/MakerMacros.h" + +DEFINE_FWK_MODULE(OverlapTagger); diff --git a/trunk/Alignment/TrackerAlignment/plugins/OverlapTagger.h b/trunk/Alignment/TrackerAlignment/plugins/OverlapTagger.h new file mode 100644 index 0000000000000..f6595a79b0674 --- /dev/null +++ b/trunk/Alignment/TrackerAlignment/plugins/OverlapTagger.h @@ -0,0 +1,51 @@ +#include <Riostream.h> + +#include "FWCore/Framework/interface/EDProducer.h" +#include "FWCore/Framework/interface/EventPrincipal.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/EventSetup.h" +#include "FWCore/Framework/interface/ESHandle.h" + +#include "FWCore/MessageLogger/interface/MessageLogger.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/ParameterSet/interface/InputTag.h" + +#include "DataFormats/Common/interface/View.h" +#include "DataFormats/SiPixelDetId/interface/PXBDetId.h" +#include "DataFormats/SiPixelDetId/interface/PXBDetId.h" +#include "DataFormats/SiPixelDetId/interface/PXFDetId.h" +#include "DataFormats/SiStripDetId/interface/TIBDetId.h" +#include "DataFormats/SiStripDetId/interface/TIDDetId.h" +#include "DataFormats/SiStripDetId/interface/TOBDetId.h" +#include "DataFormats/SiStripDetId/interface/TECDetId.h" + +#include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h" +#include "DataFormats/TrackReco/interface/Track.h" +#include "TrackingTools/PatternTools/interface/Trajectory.h" +#include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h" + +#include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2D.h" +#include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHit.h" +#include "RecoTracker/TransientTrackingRecHit/interface/TSiStripRecHit2DLocalPos.h" +#include "RecoTracker/TransientTrackingRecHit/interface/TSiPixelRecHit.h" +#include "Utilities/General/interface/ClassName.h" + +#include "DataFormats/Alignment/interface/AlignmentClusterFlag.h" +#include "DataFormats/Alignment/interface/AliClusterValueMap.h" +//#include <boost/regex.hpp> + +class OverlapTagger : public edm::EDProducer { + public: + OverlapTagger(const edm::ParameterSet &iConfig); + ~OverlapTagger(); + void produce(edm::Event &iEvent, const edm::EventSetup &iSetup); + + private: + edm::InputTag src_; + edm::InputTag srcClust_; + bool rejectBadMods_; + std::vector<unsigned int> BadModsList_; + + + int layerFromId (const DetId& id) const; +}; diff --git a/trunk/Alignment/TrackerAlignment/plugins/TreeMerger.cc b/trunk/Alignment/TrackerAlignment/plugins/TreeMerger.cc new file mode 100644 index 0000000000000..06dc15f4f6ca7 --- /dev/null +++ b/trunk/Alignment/TrackerAlignment/plugins/TreeMerger.cc @@ -0,0 +1,260 @@ +#include "Alignment/TrackerAlignment/plugins/TreeMerger.h" + +TreeMerger::TreeMerger(const edm::ParameterSet &iConfig) : + filelist_(iConfig.getParameter<string>("FileList")), + treename_(iConfig.getParameter<string>("TreeName")), + outfilename_(iConfig.getParameter<string>("OutputFile")), + // maxhits_(iConfig.getParameter<vint>("NhitsMaxLimit")) + maxhits_(iConfig.getParameter<int32_t>("NhitsMaxLimit")), + maxhitsSet_(iConfig.getParameter<edm::ParameterSet>("NhitsMaxSet")) +{ + maxPXBhits_=maxhitsSet_.getParameter<int32_t>("PXBmaxhits"); + maxPXFhits_=maxhitsSet_.getParameter<int32_t>("PXFmaxhits"); + maxTIBhits_=maxhitsSet_.getParameter<int32_t>("TIBmaxhits"); + maxTIDhits_=maxhitsSet_.getParameter<int32_t>("TIDmaxhits"); + maxTOBhits_=maxhitsSet_.getParameter<int32_t>("TOBmaxhits"); + maxTEChits_=maxhitsSet_.getParameter<int32_t>("TECmaxhits"); + //anything you want to do for initializing + cout<<"\n\n*** MAX N HITS = "<<maxhits_<<endl<<endl; +} + + +TreeMerger::~TreeMerger(){ + //default destructor + // delete out_; + // delete firsttree_; + + delete ch_; + cout<<"finished."<<endl; +} + +// ------------ method called before analyzing the first event ------------ +void TreeMerger::beginJob( const edm::EventSetup &iSetup){ + + myclock.Start(); + + //prepare the chain + ch_=new TChain(treename_.c_str()); + cout<<"The chain contains "<<ch_->GetNtrees()<<" trees"<<endl; + + //load the trees into the chain + ifstream flist(filelist_.c_str(),ios::in); + std::string filename; + std::string firstfilename; + bool first=true; + while(!flist.eof()){ + filename=""; + flist>>filename; + if(filename.empty())continue; + //cout<<"Adding "<<filename<<endl; + ch_->Add(filename.c_str()); + if(first){ + firstfilename_=filename; + first=false; + } + + } + cout<<"Now the chain contains "<<ch_->GetNtrees()<<" trees ("<<ch_->GetEntries()<<" entries)"<<endl; + + + unsigned int id_ch=0; + uint32_t nhits_ch=0,noverlaps_ch=0; + ch_->SetBranchAddress("DetId", &id_ch); + ch_->SetBranchAddress("Nhits", &nhits_ch); + ch_->SetBranchAddress("Noverlaps",&noverlaps_ch); + + ch_->SetBranchStatus("SubDet",0); + ch_->SetBranchStatus("Layer",0); + ch_->SetBranchStatus("is2D",0); + ch_->SetBranchStatus("isStereo",0); + ch_->SetBranchStatus("posX",0); + ch_->SetBranchStatus("posY",0); + ch_->SetBranchStatus("posZ",0); + ch_->SetBranchStatus("posR",0); + ch_->SetBranchStatus("posEta",0); + ch_->SetBranchStatus("posPhi",0);//now only id, nhits and noverlaps are on... + + + int totnhits(0),totnoverlaps(0); + + + //look if you find this detid in the map + DetHitMap::iterator mapiter; + DetHitMap::iterator overlapiter; + + for(int ent=0;ent<ch_->GetEntries();++ent){ + // for(int ent=0;ent<100;++ent){ + ch_->GetEntry(ent); + totnhits+=nhits_ch; + totnoverlaps+=noverlaps_ch; + + mapiter= hitmap_.find(id_ch); + if(mapiter!=hitmap_.end() ){//present, increase its value + hitmap_[id_ch]=hitmap_[id_ch]+nhits_ch; + } + else{//not present, let's add this key to the map with value=1 + hitmap_.insert(pair<uint32_t, uint32_t>(id_ch, nhits_ch)); + } + //do the same for overlaps + overlapiter= overlapmap_.find(id_ch); + if(overlapiter!=overlapmap_.end() ){ + overlapmap_[id_ch]=overlapmap_[id_ch]+noverlaps_ch; + } + else{ + overlapmap_.insert(pair<uint32_t, uint32_t>(id_ch, noverlaps_ch)); + } + + }//end loop on ent - entries in the chain + + + cout<<"Nhits in the chain: "<<totnhits<<endl; + cout<<"NOverlaps in the chain: "<<totnoverlaps<<endl; + + + myclock.Stop(); + cout<<"Finished beginJob after "<<myclock.RealTime()<<" s (real time) / "<<myclock.CpuTime()<<" s (cpu time)"<<endl; + myclock.Continue(); +}//end beginJob + + +// ------------ method called to for each event ------------ +void TreeMerger::analyze(const edm::Event&, const edm::EventSetup&){ + // cout<<firsttree_->GetEntries()<<endl; +}//end analyze + +// ------------ method called after having analyzed all the events ------------ +void TreeMerger::endJob(){ + + + //address variables in the first tree and in the chain + TFile *firstfile=new TFile(firstfilename_.c_str(),"READ"); + firsttree_=(TTree*)firstfile->Get(treename_.c_str()); + cout<<"the first tree has "<<firsttree_->GetEntries() <<" entries"<<endl; + unsigned int id=0; + uint32_t nhits=0,noverlaps=0; + float posX(-99999.0),posY(-77777.0),posZ(-88888.0); + float posEta(-6666.0),posPhi(-5555.0),posR(-4444.0); + int subdet=0; + unsigned int layer=0; + // bool is2D=false,isStereo=false; + firsttree_->SetBranchAddress("DetId", &id); + firsttree_->SetBranchAddress("Nhits", &nhits); + firsttree_->SetBranchAddress("Noverlaps",&noverlaps); + firsttree_->SetBranchAddress("SubDet", &subdet); + firsttree_->SetBranchAddress("Layer", &layer); + // firsttree_->SetBranchAddress("is2D" , &is2D); + // firsttree_->SetBranchAddress("isStereo", &isStereo); + firsttree_->SetBranchAddress("posX", &posX); + firsttree_->SetBranchAddress("posY", &posY); + firsttree_->SetBranchAddress("posZ", &posZ); + firsttree_->SetBranchAddress("posR", &posR); + firsttree_->SetBranchAddress("posEta", &posEta); + firsttree_->SetBranchAddress("posPhi", &posPhi); + + + //create and book the output + + + TFile *outfile=new TFile(outfilename_.c_str(),"RECREATE"); + out_=new TTree(treename_.c_str(),"AlignmentHitMapsTOTAL"); + unsigned int id_out=0; + uint32_t nhits_out=0,noverlaps_out=0; + float posX_out(-99999.0),posY_out(-77777.0),posZ_out(-88888.0); + float posEta_out(-6666.0),posPhi_out(-5555.0),posR_out(-4444.0); + int subdet_out=0; + unsigned int layer_out=0; + bool is2D_out=false,isStereo_out=false; + float prescfact_out=1.0; + float prescfact_overlap_out=1.0; + + out_->Branch("DetId", &id_out , "DetId/i"); + out_->Branch("Nhits", &nhits_out , "Nhits/i"); + out_->Branch("Noverlaps",&noverlaps_out,"Noverlaps/i"); + out_->Branch("SubDet", &subdet_out, "SubDet/I"); + out_->Branch("Layer", &layer_out, "Layer/i"); + out_->Branch("is2D" , &is2D_out, "is2D/B"); + out_->Branch("isStereo", &isStereo_out, "isStereo/B"); + out_->Branch("posX", &posX_out, "posX/F"); + out_->Branch("posY", &posY_out, "posY/F"); + out_->Branch("posZ", &posZ_out, "posZ/F"); + out_->Branch("posR", &posR_out, "posR/F"); + out_->Branch("posEta", &posEta_out, "posEta/F"); + out_->Branch("posPhi", &posPhi_out, "posPhi/F"); + out_->Branch("PrescaleFactor",&prescfact_out,"PrescaleFact/F"); + out_->Branch("PrescaleFactorOverlap",&prescfact_overlap_out,"PrescaleFactOverlap/F"); + + + //look if you find this detid in the map + DetHitMap::iterator mapiter; + + for(int mod=0;mod<firsttree_->GetEntries();mod++){ + //for(int mod=0;mod<100;++mod){ + // nhits_out=0; + // noverlaps_out=0; + + firsttree_->GetEntry(mod); + nhits_out=hitmap_[id]; + noverlaps_out=overlapmap_[id]; + // if(mod<25)cout<<"Nhits 1st tree: "<<nhits<<"\tTotal nhits chain: "<<nhits_out<<endl; + id_out=id; + subdet_out=subdet; + layer_out=layer; + posX_out=posX; + posY_out=posY; + posZ_out=posZ; + posR_out=posR; + posEta_out=posEta; + posPhi_out=posPhi; + + //calculate prescaling factor + int subdetmax=-1; + if(subdet_out==1)subdetmax=maxPXBhits_; + else if(subdet_out==2)subdetmax=maxPXFhits_; + else if(subdet_out==3)subdetmax=maxTIBhits_; + else if(subdet_out==4)subdetmax=maxTIDhits_; + else if(subdet_out==5)subdetmax=maxTOBhits_; + else if(subdet_out==6)subdetmax=maxTEChits_; + else subdetmax=-9999; + + if(maxhits_>-1){ + if(int(nhits_out)>maxhits_){ + prescfact_out=float(maxhits_)/float(nhits_out); + } + if(int(noverlaps_out)>maxhits_){ + prescfact_overlap_out=float(maxhits_)/float(noverlaps_out); + } + } + else if(subdetmax>0){//calculate different prescaling factors for each subdet + if(int(nhits_out)>subdetmax){ + prescfact_out=float(subdetmax/nhits_out); + } + if(int(noverlaps_out)>subdetmax){ + prescfact_overlap_out=float(subdetmax)/float(noverlaps_out); + } + } + else{ + prescfact_out=1.0; + prescfact_overlap_out=1.0; + } + out_->Fill(); + + }//end loop on mod - first tree modules + + + myclock.Stop(); + cout<<"Finished endJob after "<<myclock.RealTime()<<" s (real time) / "<<myclock.CpuTime()<<" s (cpu time)"<<endl; + cout<<"Ending the tree merging."<<endl; + out_->Write(); + cout<<"Deleting..."<<flush; + delete firstfile; + delete outfile; + + +}//end endJob + + +// ========= MODULE DEF ============== +#include "FWCore/PluginManager/interface/ModuleDef.h" +#include "FWCore/Framework/interface/MakerMacros.h" +DEFINE_FWK_MODULE(TreeMerger); + diff --git a/trunk/Alignment/TrackerAlignment/plugins/TreeMerger.h b/trunk/Alignment/TrackerAlignment/plugins/TreeMerger.h new file mode 100644 index 0000000000000..447d57ebc9c51 --- /dev/null +++ b/trunk/Alignment/TrackerAlignment/plugins/TreeMerger.h @@ -0,0 +1,59 @@ +#ifndef TrackerAlignment_TreeMerger_H +#define TrackerAlignment_TreeMerger_H + + +#include <Riostream.h> +#include <string> +#include <fstream> +#include <map> + +#include "FWCore/Framework/interface/EDAnalyzer.h" +#include "FWCore/Framework/interface/EventPrincipal.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/EventSetup.h" +#include "FWCore/Framework/interface/ESHandle.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/ParameterSet/interface/InputTag.h" + + +#include "TFile.h" +#include "TString.h" +#include "TChain.h" +#include "TStopwatch.h" + +class TreeMerger : public edm::EDAnalyzer{ + + public: + TreeMerger(const edm::ParameterSet &iConfig); + ~TreeMerger(); + void beginJob( const edm::EventSetup &iSetup); + void endJob(); + void analyze(const edm::Event&, const edm::EventSetup&); + + private: + TTree *out_;//TTree containing the merged result + TTree *firsttree_;//first tree of the list; this gives the structure to all the others + TChain *ch_;//chain containing all the tree you want to merge + std::string filelist_;//text file containing the list of input files + std::string firstfilename_; + std::string treename_;//name of the tree you want to merge (contained in the file) + std::string outfilename_;//name of the file where you want to save the output + + //Hit Population + typedef map<uint32_t,uint32_t>DetHitMap; + DetHitMap hitmap_; + DetHitMap overlapmap_; + int maxhits_;//above this number, the hit population is prescaled. Configurable for each subdet + edm::ParameterSet maxhitsSet_; + int maxPXBhits_, maxPXFhits_, maxTIBhits_, maxTIDhits_, maxTOBhits_, maxTEChits_; + + + TStopwatch myclock; + +}; + + + + +#endif diff --git a/trunk/Alignment/TrackerAlignment/python/AlignmentPrescaler_cff.py b/trunk/Alignment/TrackerAlignment/python/AlignmentPrescaler_cff.py new file mode 100644 index 0000000000000..a0a2d2dceec00 --- /dev/null +++ b/trunk/Alignment/TrackerAlignment/python/AlignmentPrescaler_cff.py @@ -0,0 +1,8 @@ +import FWCore.ParameterSet.Config as cms + +AlignmentPrescaler = cms.EDProducer("AlignmentPrescaler", + src = cms.InputTag('generalTracks'), + assomap=cms.InputTag('OverlapAssoMap'), + PrescFileName=cms.string('PrescaleFactors.root'), + PrescTreeName=cms.string('AlignmentHitMap')#if you change this be sure to be consistent with the rest of your code + ) diff --git a/trunk/Alignment/TrackerAlignment/python/OverlapTagger_cff.py b/trunk/Alignment/TrackerAlignment/python/OverlapTagger_cff.py new file mode 100644 index 0000000000000..71ace7b2f8f61 --- /dev/null +++ b/trunk/Alignment/TrackerAlignment/python/OverlapTagger_cff.py @@ -0,0 +1,9 @@ +import FWCore.ParameterSet.Config as cms + +OverlapTagger = cms.EDProducer("OverlapTagger", + src = cms.InputTag("generalTracks"), + Clustersrc = cms.InputTag("ALCARECOTkAlCosmicsCTF0T"), + rejectBadMods = cms.bool(False), + BadMods = cms.vuint32() + + )###end of module diff --git a/trunk/Alignment/TrackerAlignment/python/TreeMerger_cff.py b/trunk/Alignment/TrackerAlignment/python/TreeMerger_cff.py new file mode 100644 index 0000000000000..75691077701b8 --- /dev/null +++ b/trunk/Alignment/TrackerAlignment/python/TreeMerger_cff.py @@ -0,0 +1,16 @@ +import FWCore.ParameterSet.Config as cms + +AlignmentTreeMerger = cms.EDAnalyzer("TreeMerger", + FileList= cms.string("DQMHitMapsList.txt"), + TreeName= cms.string("AlignmentHitMaps"),#if you change this be sure to be consistent with the rest of your code + OutputFile=cms.string("AlignmentHitMapsMerged.root"), + NhitsMaxLimit=cms.int32(0),#this applies to ALL TK at the same time; no upper limit by default + NhitsMaxSet=cms.PSet(#in this way you can set different thresholds for each subdet; it is ignored if NhitsMaxLimit is higher than -1 + PXBmaxhits=cms.int32(-1),#no upper limit by default + PXFmaxhits=cms.int32(-1), + TIBmaxhits=cms.int32(-1), + TIDmaxhits=cms.int32(-1), + TOBmaxhits=cms.int32(-1), + TECmaxhits=cms.int32(-1) + ) + )