Gem cosmic muon

RecoMuon/CosmicMuonProducer/src/GEMCosmicMuon.cc

@@ -35,6 +35,39 @@
 #include "Geometry/Records/interface/MuonGeometryRecord.h"
 
 using namespace std;
+//CAS ADD MuonRecHitContainerLayered
+class MuonRecHitContainerLayered : public MuonTransientTrackingRecHit::MuonRecHitContainer {
+ protected:
+
+ public:
+  std::vector<int> Layers;
+  MuonRecHitContainerLayered( MuonTransientTrackingRecHit::MuonRecHitContainer &rechit, std::vector<int> _Layers );
+  ~MuonRecHitContainerLayered( void );
+  int GetLayer( unsigned int _L );
+  // MuonTransientTrackingRecHit::MuonRecHitContaine& GetmuRecHits(void);
+};
+
+MuonRecHitContainerLayered::MuonRecHitContainerLayered( MuonTransientTrackingRecHit::MuonRecHitContainer &rechit, std::vector<int> _Layers )
+  :MuonRecHitContainer( rechit )//MuonRecHitContaine constructor call
+{
+  Layers = _Layers;
+}
+
+MuonRecHitContainerLayered::~MuonRecHitContainerLayered( void ){
+
+}
+int MuonRecHitContainerLayered::GetLayer( unsigned int _L ){
+  try {
+    if( _L > Layers.capacity() )
+      {
+	throw  "Max number of layer exceeded";
+      }
+  }catch(const char* msg)
+    {
+      std::cerr << msg << endl;
+    }
+  return Layers[_L];  
+}
 
 class GEMCosmicMuon : public edm::stream::EDProducer<> {
 public:
@@ -53,7 +86,9 @@ class GEMCosmicMuon : public edm::stream::EDProducer<> {
   CosmicMuonSmoother* theSmoother;
   MuonServiceProxy* theService;
   KFUpdator* theUpdator;
-  unique_ptr<std::vector<TrajectorySeed> > findSeeds(MuonTransientTrackingRecHit::MuonRecHitContainer &muRecHits);
+  std::vector<int> refChambers;
+  //CAS MuonRecHitContainerLayered &muRecHits
+  unique_ptr<std::vector<TrajectorySeed> > findSeeds(MuonRecHitContainerLayered &muRecHits);
   Trajectory makeTrajectory(TrajectorySeed seed, MuonTransientTrackingRecHit::MuonRecHitContainer &muRecHits, vector<const GEMChamber*> gemChambers);
 };
 
@@ -70,6 +105,8 @@ GEMCosmicMuon::GEMCosmicMuon(const edm::ParameterSet& ps) : iev(0) {
   edm::ParameterSet smootherPSet = ps.getParameter<edm::ParameterSet>("MuonSmootherParameters");
   theSmoother = new CosmicMuonSmoother(smootherPSet,theService);
   theUpdator = new KFUpdator();
+  //CAS Add refChambers 
+  refChambers = ps.getParameter<std::vector<int>>("refChambers");//Reference chambers
   produces<reco::TrackCollection>();
   produces<TrackingRecHitCollection>();
   produces<reco::TrackExtraCollection>();
@@ -146,10 +183,15 @@ void GEMCosmicMuon::produce(edm::Event& ev, const edm::EventSetup& setup) {
     }
   }
   // cout<< "Number of GEM rechits   = " << muRecHits.size()<<endl;
-
-  trajectorySeeds =findSeeds(muRecHits);
+  /////////////////////////////////////////////////////////////////////////////////////
+  //Changed to use the new class CAS
+  /////////////////////////////////////////////////////////////////////////////////////
+  //std::vector<unsigned> Layers {, true, false};//temporal 
+  MuonRecHitContainerLayered muRecHitsl(muRecHits,refChambers);
+  trajectorySeeds =findSeeds(muRecHitsl);
+  /////////////////////////////////////////////////////////////////////////////////////
   //  cout << "GEMCosmicMuon::trajectorySeeds->size() " << trajectorySeeds->size() << endl;
-
+  
   // need to loop over seeds, make best track and save only best track
   //TrajectorySeed seed =trajectorySeeds->at(0);
   Trajectory bestTrajectory;
@@ -237,48 +279,97 @@ void GEMCosmicMuon::produce(edm::Event& ev, const edm::EventSetup& setup) {
 
 }
 
-unique_ptr<std::vector<TrajectorySeed> > GEMCosmicMuon::findSeeds(MuonTransientTrackingRecHit::MuonRecHitContainer &muRecHits)
+unique_ptr<std::vector<TrajectorySeed> > GEMCosmicMuon::findSeeds(MuonRecHitContainerLayered &muRecHits)
 {
-  unique_ptr<std::vector<TrajectorySeed> > trajectorySeeds( new vector<TrajectorySeed>());
-  for (auto hit1 : muRecHits){
-    for (auto hit2 : muRecHits){
-      if (hit1->globalPosition().y() < hit2->globalPosition().y()){
-	LocalPoint segPos = hit1->localPosition();
-	GlobalVector segDirGV(hit2->globalPosition().x() - hit1->globalPosition().x(),
-			      (hit2->globalPosition().y() - hit1->globalPosition().y()),
-			      hit2->globalPosition().z() - hit1->globalPosition().z());
-
-	segDirGV *=10;
-	LocalVector segDir = hit1->det()->toLocal(segDirGV);
-	// cout << "GEMCosmicMuon::GlobalVector " << segDirGV << endl;
-	// cout << "GEMCosmicMuon::LocalVector  " << segDir << endl;
-
-	int charge= 1;
-	LocalTrajectoryParameters param(segPos, segDir, charge);
-
-	AlgebraicSymMatrix mat(5,0);
-	mat = hit1->parametersError().similarityT( hit1->projectionMatrix() );
-	//float p_err = 0.2;
-	//mat[0][0]= p_err;
-	LocalTrajectoryError error(asSMatrix<5>(mat));
-
-	// get first hit
-	TrajectoryStateOnSurface tsos(param, error, hit1->det()->surface(), &*theService->magneticField());
-	//cout << "GEMCosmicMuon::tsos " << tsos << endl;
-	uint32_t id = hit1->rawId();
-   	PTrajectoryStateOnDet const & seedTSOS = trajectoryStateTransform::persistentState(tsos, id);
-
-	edm::OwnVector<TrackingRecHit> seedHits;
-	seedHits.push_back(hit1->hit()->clone());
-	seedHits.push_back(hit2->hit()->clone());
+  unique_ptr<std::vector<TrajectorySeed> > tmptrajectorySeeds( new vector<TrajectorySeed>());
+  for (auto hit1 : muRecHits){/* Range-based for loop, access by value, the type of hit1 is MuonTransientTrackingRecHit::MuonRecHitContaine&*/
+
+    /////////////////////////////////////////////////////////////
+    // Requesting Only seed for reference layer
+    /////////////////////////////////////////////////////////////
+    GEMDetId hit1ID(hit1->rawId());
+    //    cout<< "layer" <<hit1ID.layer() <<endl;
+    //    cout<< "chamber" <<hit1ID.chamber() <<endl;
+    bool flaghit1 = false;
+    for (auto chambers : muRecHits.Layers){
+      if( chambers % 2 == 0  ){
+	int SchamberN = chambers - 1;
+	if (SchamberN == hit1ID.chamber()){
+	  if ( hit1ID.layer() == 1 ) {
+	    flaghit1 = true;
+	  }
+	}
+      }else{
+	if ( chambers == hit1ID.chamber()){
+	  if ( hit1ID.layer() == 2 ) {
+	    flaghit1 = true;
+	  }
+	}
+      }
+
 
-	TrajectorySeed seed(seedTSOS,seedHits,alongMomentum);
-	trajectorySeeds->push_back(seed);
+      if(flaghit1){//Condition over hit1
+	bool flaghit2 = false;
+	for (auto hit2 : muRecHits){
+	  GEMDetId hit2ID(hit2->rawId());
+	  for (auto chambers : muRecHits.Layers){	 
+	    if( chambers % 2 == 0  ){
+	      int SchamberN = chambers - 1;
+	      if (SchamberN == hit2ID.chamber()){
+		if ( hit2ID.layer() == 1 ) {
+		  flaghit2 = true;
+		}
+	      }
+	    }else{
+	      if ( chambers == hit2ID.chamber()){
+		if ( hit2ID.layer() == 2 ) {
+		  flaghit2 = true;
+		}
+	      }
+	    }
+
+	    if (flaghit2){//Condition over the hit2
+	      // 	  cout<< "hit1ID_chamber" <<hit1ID.chamber() <<endl;
+	      // 	  cout<< "hit2ID_chamber" <<hit2ID.chamber() <<endl;
+	      ////////////////////////////////////////////////////////////
+
+	      if (hit1->globalPosition().y() < hit2->globalPosition().y()){/*y direction is up-down? (cosmic), see if hit1 is Higher than the hit2*/
+		LocalPoint segPos = hit1->localPosition();/*Define a LocalPoint with the position of hit1*/
+		GlobalVector segDirGV(hit2->globalPosition().x() - hit1->globalPosition().x(),/*Define a vector called segDirGV from hit2-hit1 */
+				      (hit2->globalPosition().y() - hit1->globalPosition().y()),
+				      hit2->globalPosition().z() - hit1->globalPosition().z());
+		//cout<< " hit1->localPosition().y()" << hit1->localPosition().y()<<endl;
+		//cout<< " hit1->globalPosition().y()" << hit1->globalPosition().y()<<endl;
+		segDirGV *=10;/*multiply per 10 ????*/
+		//segDirGV *=1;
+		LocalVector segDir = hit1->det()->toLocal(segDirGV);
+
+		int charge= 1;
+		LocalTrajectoryParameters param(segPos, segDir, charge);
+
+		AlgebraicSymMatrix mat(5,0);
+		mat = hit1->parametersError().similarityT( hit1->projectionMatrix() );
+		LocalTrajectoryError error(asSMatrix<5>(mat));
+
+		TrajectoryStateOnSurface tsos(param, error, hit1->det()->surface(), &*theService->magneticField());
+		uint32_t id = hit1->rawId();
+		PTrajectoryStateOnDet const & seedTSOS = trajectoryStateTransform::persistentState(tsos, id);
+
+		edm::OwnVector<TrackingRecHit> seedHits;
+		seedHits.push_back(hit1->hit()->clone());
+		seedHits.push_back(hit2->hit()->clone());
+
+		TrajectorySeed seed(seedTSOS,seedHits,alongMomentum);
+		tmptrajectorySeeds->push_back(seed);
+	      }
+	    }
+	  }
+	}
       }
     }
   }
 
-  return trajectorySeeds;
+  return tmptrajectorySeeds;
 }
 Trajectory GEMCosmicMuon::makeTrajectory(TrajectorySeed seed, MuonTransientTrackingRecHit::MuonRecHitContainer &muRecHits, vector<const GEMChamber*> gemChambers)
 {

Validation/GEMCR/BuildFile.xml

@@ -28,6 +28,7 @@
 <use name="RecoMuon/TrackingTools"/>
 <use name="RecoMuon/CosmicMuonProducer"/>
 <use name="DataFormats/GEMDigi"/>
+<use name="SimDataFormats/GEMDigiSimLink/"/>
 <export>
 <lib name="1"/>
 </export>

Validation/GEMCR/interface/gemcrValidation.h

@@ -21,6 +21,16 @@
 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
 
 
+class MuonRecHitContainerLayered : public MuonTransientTrackingRecHit::MuonRecHitContainer {
+ protected:
+
+ public:
+  std::vector<int> Layers;
+  MuonRecHitContainerLayered( MuonTransientTrackingRecHit::MuonRecHitContainer &rechit, std::vector<int> _Layers );
+  ~MuonRecHitContainerLayered( void );
+  int GetLayer( unsigned int _L );
+  // MuonTransientTrackingRecHit::MuonRecHitContaine& GetmuRecHits(void);
+};
 
 
 class gemcrValidation : public GEMBaseValidation
@@ -57,7 +67,6 @@ class gemcrValidation : public GEMBaseValidation
   std::vector<MonitorElement*> gem_chamber_bestChi2;
   std::vector<MonitorElement*> gem_chamber_track;
 
-
   MonitorElement* gemcr_g;
   MonitorElement* gem_cls_tot;
   MonitorElement* gem_bx_tot;
@@ -72,17 +81,25 @@ class gemcrValidation : public GEMBaseValidation
   MonitorElement* trajectoryh;
   MonitorElement* firedMul;
   MonitorElement* firedChamber;
-
+  MonitorElement* NumberOfSeeds;
+
 
 
   std::vector<GEMChamber> gemChambers;
   int n_ch;
   MuonServiceProxy* theService;
   CosmicMuonSmoother* theSmoother;
   KFUpdator* theUpdator;
-  std::auto_ptr<std::vector<TrajectorySeed> > findSeeds(MuonTransientTrackingRecHit::MuonRecHitContainer &muRecHits);
+  std::vector<int> refChambers;
+
+  std::auto_ptr<std::vector<TrajectorySeed> > findSeeds(MuonRecHitContainerLayered &muRecHits);
   Trajectory makeTrajectory(TrajectorySeed seed, MuonTransientTrackingRecHit::MuonRecHitContainer &muRecHits, std::vector<GEMChamber> gemChambers, GEMChamber testChamber);
   edm::EDGetToken InputTagToken_, InputTagToken_RH, InputTagToken_TR, InputTagToken_TS, InputTagToken_DG;
+
+  // ----------member data ---------------------------      
+
+  //   static unsigned int SeedsNumber; 
+
 };
 
 #endif