CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_6_1_1/src/FastSimulation/Tracking/plugins/TrajectorySeedProducer.cc

Go to the documentation of this file.
00001 #include <memory>
00002 
00003 #include "FWCore/Framework/interface/Event.h"
00004 #include "FWCore/Framework/interface/EventSetup.h"
00005 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00006 #include "FWCore/Framework/interface/ESHandle.h"
00007 
00008 #include "DataFormats/Common/interface/Handle.h"
00009 #include "DataFormats/Common/interface/OwnVector.h"
00010 #include "DataFormats/TrackerRecHit2D/interface/SiTrackerGSRecHit2DCollection.h" 
00011 #include "DataFormats/TrackerRecHit2D/interface/SiTrackerGSMatchedRecHit2DCollection.h" 
00012 #include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h"
00013 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
00014 #include "DataFormats/VertexReco/interface/Vertex.h"
00015 
00016 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
00017 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
00018 
00019 #include "FastSimulation/ParticlePropagator/interface/MagneticFieldMapRecord.h"
00020 
00021 #include "FastSimulation/Tracking/plugins/TrajectorySeedProducer.h"
00022 #include "FastSimulation/Tracking/interface/TrackerRecHit.h"
00023 
00024 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
00025 #include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
00026 
00027 #include "TrackingTools/TrajectoryParametrization/interface/CurvilinearTrajectoryError.h"
00028 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
00029 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
00030 
00031 #include "Geometry/CommonDetUnit/interface/GeomDetUnit.h"
00032 #include "DataFormats/DetId/interface/DetId.h"
00033 
00034 #include "FastSimulation/BaseParticlePropagator/interface/BaseParticlePropagator.h"
00035 #include "FastSimulation/ParticlePropagator/interface/ParticlePropagator.h"
00036 
00037 //Propagator withMaterial
00038 #include "TrackingTools/MaterialEffects/interface/PropagatorWithMaterial.h"
00039 //analyticalpropagator
00040 //#include "TrackingTools/GeomPropagators/interface/AnalyticalPropagator.h"
00041 
00042 
00043 //
00044 
00045 //for debug only 
00046 //#define FAMOS_DEBUG
00047 
00048 TrajectorySeedProducer::TrajectorySeedProducer(const edm::ParameterSet& conf) :thePropagator(0)
00049 {  
00050 
00051   // The input tag for the beam spot
00052   theBeamSpot = conf.getParameter<edm::InputTag>("beamSpot");
00053 
00054   // The name of the TrajectorySeed Collections
00055   seedingAlgo = conf.getParameter<std::vector<std::string> >("seedingAlgo");
00056   for ( unsigned i=0; i<seedingAlgo.size(); ++i )
00057     produces<TrajectorySeedCollection>(seedingAlgo[i]);
00058 
00059   // The smallest true pT for a track candidate
00060   pTMin = conf.getParameter<std::vector<double> >("pTMin");
00061   if ( pTMin.size() != seedingAlgo.size() ) 
00062     throw cms::Exception("FastSimulation/TrajectorySeedProducer ") 
00063       << " WARNING : pTMin does not have the proper size "
00064       << std::endl;
00065 
00066   for ( unsigned i=0; i<pTMin.size(); ++i )
00067     pTMin[i] *= pTMin[i];  // Cut is done of perp2() - CPU saver
00068   
00069   // The smallest number of Rec Hits for a track candidate
00070   minRecHits = conf.getParameter<std::vector<unsigned int> >("minRecHits");
00071   if ( minRecHits.size() != seedingAlgo.size() ) 
00072     throw cms::Exception("FastSimulation/TrajectorySeedProducer ") 
00073       << " WARNING : minRecHits does not have the proper size "
00074       << std::endl;
00075   // Set the overall number hits to be checked
00076   absMinRecHits = 0;
00077   for ( unsigned ialgo=0; ialgo<minRecHits.size(); ++ialgo ) 
00078     if ( minRecHits[ialgo] > absMinRecHits ) absMinRecHits = minRecHits[ialgo];
00079 
00080   // The smallest true impact parameters (d0 and z0) for a track candidate
00081   maxD0 = conf.getParameter<std::vector<double> >("maxD0");
00082   if ( maxD0.size() != seedingAlgo.size() ) 
00083     throw cms::Exception("FastSimulation/TrajectorySeedProducer ") 
00084       << " WARNING : maxD0 does not have the proper size "
00085       << std::endl;
00086 
00087   maxZ0 = conf.getParameter<std::vector<double> >("maxZ0");
00088   if ( maxZ0.size() != seedingAlgo.size() ) 
00089     throw cms::Exception("FastSimulation/TrajectorySeedProducer ") 
00090       << " WARNING : maxZ0 does not have the proper size "
00091       << std::endl;
00092 
00093   // The name of the hit producer
00094   hitProducer = conf.getParameter<edm::InputTag>("HitProducer");
00095 
00096   // The cuts for seed cleaning
00097   seedCleaning = conf.getParameter<bool>("seedCleaning");
00098 
00099   // Number of hits needed for a seed
00100   numberOfHits = conf.getParameter<std::vector<unsigned int> >("numberOfHits");
00101   if ( numberOfHits.size() != seedingAlgo.size() ) 
00102     throw cms::Exception("FastSimulation/TrajectorySeedProducer ") 
00103       << " WARNING : numberOfHits does not have the proper size "
00104       << std::endl;
00105 
00106   // Seeding based on muons
00107   selectMuons = conf.getParameter<bool>("selectMuons");
00108 
00109   // Layers
00110   newSyntax = conf.getParameter<bool>("newSyntax");
00111   if (newSyntax) {
00112     layerList = conf.getParameter<std::vector<std::string> >("layerList");
00113     // (AG)  for (unsigned i=0; i<layerList.size();i++) std::cout << "------- Layers = " << layerList[i] << std::endl;
00114   } else {
00115     // TO BE DELETED (AG)
00116     firstHitSubDetectorNumber = 
00117       conf.getParameter<std::vector<unsigned int> >("firstHitSubDetectorNumber");
00118     if ( firstHitSubDetectorNumber.size() != seedingAlgo.size() ) 
00119       throw cms::Exception("FastSimulation/TrajectorySeedProducer ") 
00120         << " WARNING : firstHitSubDetectorNumber does not have the proper size "
00121         << std::endl;
00122     
00123     std::vector<unsigned int> firstSubDets = 
00124       conf.getParameter<std::vector<unsigned int> >("firstHitSubDetectors");
00125     unsigned isub1 = 0;
00126     unsigned check1 = 0;
00127     firstHitSubDetectors.resize(seedingAlgo.size());
00128     for ( unsigned ialgo=0; ialgo<firstHitSubDetectorNumber.size(); ++ialgo ) { 
00129       check1 += firstHitSubDetectorNumber[ialgo];
00130       for ( unsigned idet=0; idet<firstHitSubDetectorNumber[ialgo]; ++idet ) { 
00131         firstHitSubDetectors[ialgo].push_back(firstSubDets[isub1++]);
00132       }
00133     }
00134     if ( firstSubDets.size() != check1 ) 
00135       throw cms::Exception("FastSimulation/TrajectorySeedProducer ") 
00136         << " WARNING : firstHitSubDetectors does not have the proper size (should be " << check1 << ")"
00137         << std::endl;
00138     
00139     
00140     secondHitSubDetectorNumber = 
00141       conf.getParameter<std::vector<unsigned int> >("secondHitSubDetectorNumber");
00142     if ( secondHitSubDetectorNumber.size() != seedingAlgo.size() ) 
00143       throw cms::Exception("FastSimulation/TrajectorySeedProducer ") 
00144         << " WARNING : secondHitSubDetectorNumber does not have the proper size "
00145         << std::endl;
00146     
00147     std::vector<unsigned int> secondSubDets = 
00148       conf.getParameter<std::vector<unsigned int> >("secondHitSubDetectors");
00149     unsigned isub2 = 0;
00150     unsigned check2 = 0;
00151     secondHitSubDetectors.resize(seedingAlgo.size());
00152     for ( unsigned ialgo=0; ialgo<secondHitSubDetectorNumber.size(); ++ialgo ) { 
00153       check2 += secondHitSubDetectorNumber[ialgo];
00154       for ( unsigned idet=0; idet<secondHitSubDetectorNumber[ialgo]; ++idet ) { 
00155         secondHitSubDetectors[ialgo].push_back(secondSubDets[isub2++]);
00156       }
00157     }
00158     if ( secondSubDets.size() != check2 ) 
00159       throw cms::Exception("FastSimulation/TrajectorySeedProducer ") 
00160         << " WARNING : secondHitSubDetectors does not have the proper size (should be " << check2 << ")"
00161         << std::endl;
00162     
00163     thirdHitSubDetectorNumber = 
00164       conf.getParameter<std::vector<unsigned int> >("thirdHitSubDetectorNumber");
00165     if ( thirdHitSubDetectorNumber.size() != seedingAlgo.size() ) 
00166       throw cms::Exception("FastSimulation/TrajectorySeedProducer ") 
00167         << " WARNING : thirdHitSubDetectorNumber does not have the proper size "
00168         << std::endl;
00169     
00170     std::vector<unsigned int> thirdSubDets = 
00171       conf.getParameter<std::vector<unsigned int> >("thirdHitSubDetectors");
00172     unsigned isub3 = 0;
00173     unsigned check3 = 0;
00174     thirdHitSubDetectors.resize(seedingAlgo.size());
00175     for ( unsigned ialgo=0; ialgo<thirdHitSubDetectorNumber.size(); ++ialgo ) { 
00176       check3 += thirdHitSubDetectorNumber[ialgo];
00177       for ( unsigned idet=0; idet<thirdHitSubDetectorNumber[ialgo]; ++idet ) { 
00178         thirdHitSubDetectors[ialgo].push_back(thirdSubDets[isub3++]);
00179       }
00180     }
00181     if ( thirdSubDets.size() != check3 ) 
00182       throw cms::Exception("FastSimulation/TrajectorySeedProducer ") 
00183         << " WARNING : thirdHitSubDetectors does not have the proper size (should be " << check3 << ")"
00184         << std::endl;
00185     
00186     originRadius = conf.getParameter<std::vector<double> >("originRadius");
00187     if ( originRadius.size() != seedingAlgo.size() ) 
00188       throw cms::Exception("FastSimulation/TrajectorySeedProducer ") 
00189         << " WARNING : originRadius does not have the proper size "
00190         << std::endl;
00191     
00192     originHalfLength = conf.getParameter<std::vector<double> >("originHalfLength");
00193     if ( originHalfLength.size() != seedingAlgo.size() ) 
00194       throw cms::Exception("FastSimulation/TrajectorySeedProducer ") 
00195         << " WARNING : originHalfLength does not have the proper size "
00196         << std::endl;
00197     
00198     originpTMin = conf.getParameter<std::vector<double> >("originpTMin");
00199     if ( originpTMin.size() != seedingAlgo.size() ) 
00200       throw cms::Exception("FastSimulation/TrajectorySeedProducer ") 
00201         << " WARNING : originpTMin does not have the proper size "
00202         << std::endl;
00203     
00204     primaryVertices = conf.getParameter<std::vector<edm::InputTag> >("primaryVertices");
00205     if ( primaryVertices.size() != seedingAlgo.size() ) 
00206       throw cms::Exception("FastSimulation/TrajectorySeedProducer ") 
00207         << " WARNING : primaryVertices does not have the proper size "
00208         << std::endl;
00209     
00210     zVertexConstraint = conf.getParameter<std::vector<double> >("zVertexConstraint");
00211     if ( zVertexConstraint.size() != seedingAlgo.size() ) 
00212       throw cms::Exception("FastSimulation/TrajectorySeedProducer ") 
00213         << " WARNING : zVertexConstraint does not have the proper size "
00214         << std::endl;
00215   }
00216 
00217 }
00218 
00219   
00220 // Virtual destructor needed.
00221 TrajectorySeedProducer::~TrajectorySeedProducer() {
00222   
00223   if(thePropagator) delete thePropagator;
00224 
00225   // do nothing
00226 #ifdef FAMOS_DEBUG
00227   std::cout << "TrajectorySeedProducer destructed" << std::endl;
00228 #endif
00229 
00230 } 
00231  
00232 void 
00233 TrajectorySeedProducer::beginRun(edm::Run & run, const edm::EventSetup & es) {
00234 
00235   //services
00236   //  es.get<TrackerRecoGeometryRecord>().get(theGeomSearchTracker);
00237 
00238   edm::ESHandle<MagneticField>          magField;
00239   edm::ESHandle<TrackerGeometry>        geometry;
00240   edm::ESHandle<MagneticFieldMap>       magFieldMap;
00241 
00242 
00243   es.get<IdealMagneticFieldRecord>().get(magField);
00244   es.get<TrackerDigiGeometryRecord>().get(geometry);
00245   es.get<MagneticFieldMapRecord>().get(magFieldMap);
00246 
00247   theMagField = &(*magField);
00248   theGeometry = &(*geometry);
00249   theFieldMap = &(*magFieldMap);
00250 
00251   thePropagator = new PropagatorWithMaterial(alongMomentum,0.105,&(*theMagField)); 
00252 
00253   const GlobalPoint g(0.,0.,0.);
00254 
00255 }
00256   
00257   // Functions that gets called by framework every event
00258 void 
00259 TrajectorySeedProducer::produce(edm::Event& e, const edm::EventSetup& es) {        
00260 
00261 
00262   //  if( seedingAlgo[0] ==  "FourthPixelLessPairs") std::cout << "Seed producer in 4th iteration " << std::endl;
00263 
00264 #ifdef FAMOS_DEBUG
00265   std::cout << "################################################################" << std::endl;
00266   std::cout << " TrajectorySeedProducer produce init " << std::endl;
00267 #endif
00268 
00269 
00270 
00271   unsigned nSimTracks = 0;
00272   unsigned nTracksWithHits = 0;
00273   unsigned nTracksWithPT = 0;
00274   unsigned nTracksWithD0Z0 = 0;
00275   //  unsigned nTrackCandidates = 0;
00276   PTrajectoryStateOnDet initialState;
00277   
00278   // Output
00279   std::vector<TrajectorySeedCollection*>
00280     output(seedingAlgo.size(),static_cast<TrajectorySeedCollection*>(0));
00281   for ( unsigned ialgo=0; ialgo<seedingAlgo.size(); ++ialgo ) { 
00282     //    std::auto_ptr<TrajectorySeedCollection> p(new TrajectorySeedCollection );
00283     output[ialgo] = new TrajectorySeedCollection;
00284   }
00285   
00286   // Beam spot
00287   edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
00288   e.getByLabel(theBeamSpot,recoBeamSpotHandle); 
00289   math::XYZPoint BSPosition_ = recoBeamSpotHandle->position();
00290 
00291   //not used anymore. take the value from the py
00292 
00293   //double sigmaZ=recoBeamSpotHandle->sigmaZ();
00294   //double sigmaZ0Error=recoBeamSpotHandle->sigmaZ0Error();
00295   //double sigmaz0=sqrt(sigmaZ*sigmaZ+sigmaZ0Error*sigmaZ0Error);
00296   x0 = BSPosition_.X();
00297   y0 = BSPosition_.Y();
00298   z0 = BSPosition_.Z();
00299 
00300   // SimTracks and SimVertices
00301   edm::Handle<edm::SimTrackContainer> theSimTracks;
00302   e.getByLabel("famosSimHits",theSimTracks);
00303   
00304   edm::Handle<edm::SimVertexContainer> theSimVtx;
00305   e.getByLabel("famosSimHits",theSimVtx);
00306 
00307 #ifdef FAMOS_DEBUG
00308   std::cout << " Step A: SimTracks found " << theSimTracks->size() << std::endl;
00309 #endif
00310   
00311   //  edm::Handle<SiTrackerGSRecHit2DCollection> theGSRecHits;
00312   edm::Handle<SiTrackerGSMatchedRecHit2DCollection> theGSRecHits;
00313   e.getByLabel(hitProducer, theGSRecHits);
00314   
00315   // No tracking attempted if no hits (but put an empty collection in the event)!
00316 #ifdef FAMOS_DEBUG
00317   std::cout << " Step B: Full GS RecHits found " << theGSRecHits->size() << std::endl;
00318 #endif
00319   if(theGSRecHits->size() == 0) {
00320     for ( unsigned ialgo=0; ialgo<seedingAlgo.size(); ++ialgo ) {
00321       std::auto_ptr<TrajectorySeedCollection> p(output[ialgo]);
00322       e.put(p,seedingAlgo[ialgo]);
00323     }
00324     return;
00325   }
00326   
00327   // Primary vertices
00328   vertices = std::vector<const reco::VertexCollection*>
00329     (seedingAlgo.size(),static_cast<const reco::VertexCollection*>(0));
00330   for ( unsigned ialgo=0; ialgo<seedingAlgo.size(); ++ialgo ) { 
00331     //PAT Attempt!!!! 
00332 
00333    //originHalfLength[ialgo] = 3.*sigmaz0; // Overrides the configuration
00334     edm::Handle<reco::VertexCollection> aHandle;
00335     bool isVertexCollection = e.getByLabel(primaryVertices[ialgo],aHandle);
00336     if (!isVertexCollection ) continue;
00337     vertices[ialgo] = &(*aHandle);
00338   }
00339   
00340 #ifdef FAMOS_DEBUG
00341   std::cout << " Step C: Loop over the RecHits, track by track " << std::endl;
00342 #endif
00343 
00344   // The vector of simTrack Id's carrying GSRecHits
00345   const std::vector<unsigned> theSimTrackIds = theGSRecHits->ids();
00346 
00347   // loop over SimTrack Id's
00348   for ( unsigned tkId=0;  tkId != theSimTrackIds.size(); ++tkId ) {
00349 
00350 #ifdef FAMOS_DEBUG
00351     std::cout << "Track number " << tkId << std::endl;
00352 #endif
00353 
00354     ++nSimTracks;
00355     unsigned simTrackId = theSimTrackIds[tkId];
00356     const SimTrack& theSimTrack = (*theSimTracks)[simTrackId]; 
00357 #ifdef FAMOS_DEBUG
00358     std::cout << "Pt = " << std::sqrt(theSimTrack.momentum().Perp2()) 
00359               << " eta " << theSimTrack.momentum().Eta()
00360               << " pdg ID " << theSimTrack.type()
00361               << std::endl;
00362 #endif
00363 
00364     // Select only muons, if requested
00365     if (selectMuons && abs(theSimTrack.type()) != 13) continue;
00366     
00367     // Check that the sim track comes from the main vertex (loose cut)
00368     int vertexIndex = theSimTrack.vertIndex();
00369     const SimVertex& theSimVertex = (*theSimVtx)[vertexIndex]; 
00370 #ifdef FAMOS_DEBUG
00371     std::cout << " o SimTrack " << theSimTrack << std::endl;
00372     std::cout << " o SimVertex " << theSimVertex << std::endl;
00373 #endif
00374     
00375     BaseParticlePropagator theParticle = 
00376       BaseParticlePropagator( 
00377          RawParticle(XYZTLorentzVector(theSimTrack.momentum().px(),
00378                                        theSimTrack.momentum().py(),
00379                                        theSimTrack.momentum().pz(),
00380                                        theSimTrack.momentum().e()),
00381                      XYZTLorentzVector(theSimVertex.position().x(),
00382                                        theSimVertex.position().y(),
00383                                        theSimVertex.position().z(),
00384                                        theSimVertex.position().t())),
00385                      0.,0.,4.);
00386     theParticle.setCharge((*theSimTracks)[simTrackId].charge());
00387 
00388     SiTrackerGSMatchedRecHit2DCollection::range theRecHitRange = theGSRecHits->get(simTrackId);
00389     SiTrackerGSMatchedRecHit2DCollection::const_iterator theRecHitRangeIteratorBegin = theRecHitRange.first;
00390     SiTrackerGSMatchedRecHit2DCollection::const_iterator theRecHitRangeIteratorEnd   = theRecHitRange.second;
00391     SiTrackerGSMatchedRecHit2DCollection::const_iterator iterRecHit;
00392     SiTrackerGSMatchedRecHit2DCollection::const_iterator iterRecHit1;
00393     SiTrackerGSMatchedRecHit2DCollection::const_iterator iterRecHit2;
00394     SiTrackerGSMatchedRecHit2DCollection::const_iterator iterRecHit3;
00395 
00396     // Check the number of layers crossed
00397     unsigned numberOfRecHits = 0;
00398     TrackerRecHit previousHit, currentHit;
00399     for ( iterRecHit = theRecHitRangeIteratorBegin; 
00400           iterRecHit != theRecHitRangeIteratorEnd; 
00401           ++iterRecHit) { 
00402       previousHit = currentHit;
00403       currentHit = TrackerRecHit(&(*iterRecHit),theGeometry);
00404       if ( currentHit.isOnTheSameLayer(previousHit) ) continue;
00405       ++numberOfRecHits;
00406       if ( numberOfRecHits == absMinRecHits ) break;
00407     }
00408 
00409     // Loop on the successive seedings
00410     for ( unsigned int ialgo = 0; ialgo < seedingAlgo.size(); ++ialgo ) { 
00411 
00412 #ifdef FAMOS_DEBUG
00413       std::cout << "Algo " << seedingAlgo[ialgo] << std::endl;
00414 #endif
00415 
00416       // Request a minimum number of RecHits for the track to give a seed.
00417 #ifdef FAMOS_DEBUG
00418       std::cout << "The number of RecHits = " << numberOfRecHits << std::endl;
00419 #endif
00420       if ( numberOfRecHits < minRecHits[ialgo] ) continue;
00421       ++nTracksWithHits;
00422 
00423       // Request a minimum pT for the sim track
00424       if ( theSimTrack.momentum().Perp2() < pTMin[ialgo] ) continue;
00425       ++nTracksWithPT;
00426       
00427       // Cut on sim track impact parameters
00428       if ( theParticle.xyImpactParameter(x0,y0) > maxD0[ialgo] ) continue;
00429       if ( fabs( theParticle.zImpactParameter(x0,y0) - z0 ) > maxZ0[ialgo] ) continue;
00430       ++nTracksWithD0Z0;
00431       
00432       std::vector<TrackerRecHit > 
00433         theSeedHits(numberOfHits[ialgo],
00434                     static_cast<TrackerRecHit >(TrackerRecHit()));
00435       TrackerRecHit& theSeedHits0 = theSeedHits[0];
00436       TrackerRecHit& theSeedHits1 = theSeedHits[1];
00437       TrackerRecHit& theSeedHits2 = theSeedHits[2];
00438       bool compatible = false;
00439       for ( iterRecHit1 = theRecHitRangeIteratorBegin; iterRecHit1 != theRecHitRangeIteratorEnd; ++iterRecHit1) {
00440         theSeedHits[0] = TrackerRecHit(&(*iterRecHit1),theGeometry);
00441 #ifdef FAMOS_DEBUG
00442         std::cout << "The first hit position = " << theSeedHits0.globalPosition() << std::endl;
00443         std::cout << "The first hit subDetId = " << theSeedHits0.subDetId() << std::endl;
00444         std::cout << "The first hit layer    = " << theSeedHits0.layerNumber() << std::endl;
00445 #endif
00446 
00447         // Check if inside the requested detectors
00448         bool isInside = true;
00449         if (!selectMuons) {
00450           if (newSyntax) 
00451             isInside = true; // AG placeholder
00452           else
00453             isInside = theSeedHits0.subDetId() < firstHitSubDetectors[ialgo][0];
00454           //    bool isInside = theSeedHits0.subDetId() < firstHitSubDetectors[ialgo][0];
00455           if ( isInside ) continue;
00456         }
00457 
00458         // Check if on requested detectors
00459         //      bool isOndet =  theSeedHits0.isOnRequestedDet(firstHitSubDetectors[ialgo]);
00460         bool isOndet = true;
00461         if (!selectMuons) {
00462           if (newSyntax) 
00463             isOndet = theSeedHits0.isOnRequestedDet(layerList);
00464           else
00465             isOndet = theSeedHits0.isOnRequestedDet(firstHitSubDetectors[ialgo], seedingAlgo[ialgo]);
00466           //    bool isOndet =  theSeedHits0.isOnRequestedDet(firstHitSubDetectors[ialgo], seedingAlgo[ialgo]);
00467           //    if ( !isOndet ) break;
00468           if ( !isOndet ) continue;
00469         }
00470 
00471 #ifdef FAMOS_DEBUG
00472         std::cout << "Apparently the first hit is on the requested detector! " << std::endl;
00473 #endif
00474         for ( iterRecHit2 = iterRecHit1+1; iterRecHit2 != theRecHitRangeIteratorEnd; ++iterRecHit2) {
00475           theSeedHits[1] = TrackerRecHit(&(*iterRecHit2),theGeometry);
00476 #ifdef FAMOS_DEBUG
00477           std::cout << "The second hit position = " << theSeedHits1.globalPosition() << std::endl;
00478           std::cout << "The second hit subDetId = " << theSeedHits1.subDetId() << std::endl;
00479           std::cout << "The second hit layer    = " << theSeedHits1.layerNumber() << std::endl;
00480 #endif
00481 
00482           if (!selectMuons) {
00483             // Check if inside the requested detectors
00484             if (newSyntax) 
00485               isInside = true; // AG placeholder
00486             else
00487               isInside = theSeedHits1.subDetId() < secondHitSubDetectors[ialgo][0];
00488             if ( isInside ) continue;
00489 
00490             // Check if on requested detectors
00491             if (newSyntax) 
00492               isOndet = theSeedHits1.isOnRequestedDet(layerList);
00493             else
00494               isOndet =  theSeedHits1.isOnRequestedDet(secondHitSubDetectors[ialgo], seedingAlgo[ialgo]);
00495             if ( !isOndet ) break;
00496           }
00497 
00498           // Check if on the same layer as previous hit
00499           if ( theSeedHits1.isOnTheSameLayer(theSeedHits0) ) continue;
00500 
00501 #ifdef FAMOS_DEBUG
00502           std::cout << "Apparently the second hit is on the requested detector! " << std::endl;
00503 #endif
00504           GlobalPoint gpos1 = theSeedHits0.globalPosition();
00505           GlobalPoint gpos2 = theSeedHits1.globalPosition();
00506           bool forward = theSeedHits0.isForward();
00507           double error = std::sqrt(theSeedHits0.largerError()+theSeedHits1.largerError());
00508           //      compatible = compatibleWithVertex(gpos1,gpos2,ialgo);
00509           //added out of desperation    
00510           if(seedingAlgo[0] == "PixelLess" ||  seedingAlgo[0] ==  "TobTecLayerPairs"){
00511             compatible = true;
00512             //std::cout << "Algo " << seedingAlgo[0] << "Det/layer = " << theSeedHits0.subDetId() << "/" <<  theSeedHits0.layerNumber() << std::endl;
00513           } else {
00514             compatible = compatibleWithBeamAxis(gpos1,gpos2,error,forward,ialgo);
00515           }
00516           
00517 #ifdef FAMOS_DEBUG
00518           std::cout << "Algo" << seedingAlgo[0] << "\t Are the two hits compatible with the PV? " << compatible << std::endl;
00519 #endif
00520 
00521           if (!selectMuons) {
00522             // Check if the pair is on the requested dets
00523             if ( numberOfHits[ialgo] == 2 ) {
00524               
00525               if ( seedingAlgo[0] ==  "ThirdMixedPairs" ){
00526                 compatible = compatible && theSeedHits[0].makesAPairWith3rd(theSeedHits[1]);
00527               } else {
00528                 compatible = compatible && theSeedHits[0].makesAPairWith(theSeedHits[1]);
00529                 //check
00530                 /*
00531                   if((seedingAlgo[0] == "PixelLess" ||  seedingAlgo[0] ==  "TobTecLayerPairs") && !compatible) 
00532                   std::cout << "NOT Compatible " <<  seedingAlgo[0] 
00533                   <<  "Hit 1 Det/layer/ring = " << theSeedHits0.subDetId() << "/" <<  theSeedHits0.layerNumber() << "/" << theSeedHits0.ringNumber() 
00534                   <<  "\tHit 2 Det/layer/ring = " << theSeedHits1.subDetId() << "/" <<  theSeedHits1.layerNumber() << "/" << theSeedHits1.ringNumber() <<  std::endl;
00535                 */
00536               }
00537             }   
00538           }    
00539           
00540           // Reject non suited pairs
00541           if ( !compatible ) continue;
00542 
00543 #ifdef FAMOS_DEBUG
00544           std::cout << "Pair kept! " << std::endl;
00545 #endif
00546 
00547           // Leave here if only two hits are required.
00548           if ( numberOfHits[ialgo] == 2 ) break; 
00549           
00550           compatible = false;
00551           // Check if there is a third satisfying hit otherwise
00552           for ( iterRecHit3 = iterRecHit2+1; iterRecHit3 != theRecHitRangeIteratorEnd; ++iterRecHit3) {
00553             theSeedHits[2] = TrackerRecHit(&(*iterRecHit3),theGeometry);
00554 #ifdef FAMOS_DEBUG
00555             std::cout << "The third hit position = " << theSeedHits2.globalPosition() << std::endl;
00556             std::cout << "The third hit subDetId = " << theSeedHits2.subDetId() << std::endl;
00557             std::cout << "The third hit layer    = " << theSeedHits2.layerNumber() << std::endl;
00558 #endif
00559 
00560             if (!selectMuons) {
00561               // Check if inside the requested detectors
00562               if (newSyntax) 
00563                 isInside = true; // AG placeholder
00564               else 
00565                 isInside = theSeedHits2.subDetId() < thirdHitSubDetectors[ialgo][0];
00566               if ( isInside ) continue;
00567             
00568               // Check if on requested detectors
00569               if (newSyntax) 
00570                 isOndet = theSeedHits2.isOnRequestedDet(layerList);
00571               else 
00572                 isOndet =  theSeedHits2.isOnRequestedDet(thirdHitSubDetectors[ialgo], seedingAlgo[ialgo]);
00573               //            if ( !isOndet ) break;
00574               if ( !isOndet ) continue;
00575             }
00576 
00577             // Check if on the same layer as previous hit
00578             compatible = !(theSeedHits2.isOnTheSameLayer(theSeedHits1));
00579 
00580             // Check if the triplet is on the requested det combination
00581             if (!selectMuons) compatible = compatible && theSeedHits[0].makesATripletWith(theSeedHits[1],theSeedHits[2]);
00582 
00583 #ifdef FAMOS_DEBUG
00584             if ( compatible ) 
00585               std::cout << "Apparently the third hit is on the requested detector! " << std::endl;
00586 #endif
00587 
00588             if ( compatible ) break;      
00589 
00590           }
00591 
00592           if ( compatible ) break;
00593 
00594         }
00595 
00596         if ( compatible ) break;
00597 
00598       }
00599 
00600       // There is no compatible seed for this track with this seeding algorithm 
00601       // Go to next algo
00602       if ( !compatible ) continue;
00603 #ifdef FAMOS_DEBUG
00604       if ( compatible ) 
00605         std::cout << "@@@ There is at least a compatible seed" << std::endl;
00606       else
00607         std::cout << "@@@ There is no compatible seed" << std::endl;
00608 #endif
00609       
00610 
00611 #ifdef FAMOS_DEBUG
00612       std::cout << "Preparing to create the TrajectorySeed" << std::endl;
00613 #endif
00614       // The seed is validated -> include in the collection
00615       // 1) Create the vector of RecHits
00616       edm::OwnVector<TrackingRecHit> recHits;
00617       for ( unsigned ih=0; ih<theSeedHits.size(); ++ih ) {
00618         TrackingRecHit* aTrackingRecHit = theSeedHits[ih].hit()->clone();
00619         recHits.push_back(aTrackingRecHit);
00620       }
00621 #ifdef FAMOS_DEBUG
00622       std::cout << "with " << recHits.size() << " hits." << std::endl;
00623 #endif
00624 
00625       // 2) Create the initial state
00626       //   a) origin vertex
00627       GlobalPoint  position((*theSimVtx)[vertexIndex].position().x(),
00628                             (*theSimVtx)[vertexIndex].position().y(),
00629                             (*theSimVtx)[vertexIndex].position().z());
00630       
00631       //   b) initial momentum
00632       GlobalVector momentum( (*theSimTracks)[simTrackId].momentum().x() , 
00633                              (*theSimTracks)[simTrackId].momentum().y() , 
00634                              (*theSimTracks)[simTrackId].momentum().z() );
00635       //   c) electric charge
00636       float        charge   = (*theSimTracks)[simTrackId].charge();
00637       //  -> inital parameters
00638       GlobalTrajectoryParameters initialParams(position,momentum,(int)charge,theMagField);
00639       //  -> large initial errors
00640       AlgebraicSymMatrix55 errorMatrix= AlgebraicMatrixID();      
00641       // errorMatrix = errorMatrix * 10;
00642 
00643       //this line help the fit succeed in the case of pixelless tracks (4th and 5th iteration)
00644       //for the future: probably the best thing is to use the mini-kalmanFilter
00645       if(theSeedHits0.subDetId() !=1 || theSeedHits0.subDetId() !=2) errorMatrix = errorMatrix * 0.0000001;
00646 
00647 
00648 
00649 #ifdef FAMOS_DEBUG
00650       std::cout << "TrajectorySeedProducer: SimTrack parameters " << std::endl;
00651       std::cout << "\t\t pT  = " << (*theSimTracks)[simTrackId].momentum().Pt() << std::endl;
00652       std::cout << "\t\t eta = " << (*theSimTracks)[simTrackId].momentum().Eta()  << std::endl;
00653       std::cout << "\t\t phi = " << (*theSimTracks)[simTrackId].momentum().Phi()  << std::endl;
00654       std::cout << "TrajectorySeedProducer: AlgebraicSymMatrix " << errorMatrix << std::endl;
00655 #endif
00656       CurvilinearTrajectoryError initialError(errorMatrix);
00657       // -> initial state
00658       FreeTrajectoryState initialFTS(initialParams, initialError);      
00659 #ifdef FAMOS_DEBUG
00660       std::cout << "TrajectorySeedProducer: FTS momentum " << initialFTS.momentum() << std::endl;
00661 #endif
00662       // const GeomDetUnit* initialLayer = theGeometry->idToDetUnit( recHits.front().geographicalId() );
00663       const GeomDet* initialLayer = theGeometry->idToDet( recHits.front().geographicalId() );
00664 
00665       //this is wrong because the FTS is defined at vertex, and it need to be properly propagated.
00666       //      const TrajectoryStateOnSurface initialTSOS(initialFTS, initialLayer->surface());      
00667 
00668       const TrajectoryStateOnSurface initialTSOS = thePropagator->propagate(initialFTS,initialLayer->surface()) ;
00669       if (!initialTSOS.isValid()) continue; 
00670 
00671 #ifdef FAMOS_DEBUG
00672       std::cout << "TrajectorySeedProducer: TSOS global momentum "    << initialTSOS.globalMomentum() << std::endl;
00673       std::cout << "\t\t\tpT = "                                     << initialTSOS.globalMomentum().perp() << std::endl;
00674       std::cout << "\t\t\teta = "                                    << initialTSOS.globalMomentum().eta() << std::endl;
00675       std::cout << "\t\t\tphi = "                                    << initialTSOS.globalMomentum().phi() << std::endl;
00676       std::cout << "TrajectorySeedProducer: TSOS local momentum "     << initialTSOS.localMomentum()  << std::endl;
00677       std::cout << "TrajectorySeedProducer: TSOS local error "        << initialTSOS.localError().positionError() << std::endl;
00678       std::cout << "TrajectorySeedProducer: TSOS local error matrix " << initialTSOS.localError().matrix() << std::endl;
00679       std::cout << "TrajectorySeedProducer: TSOS surface side "       << initialTSOS.surfaceSide()    << std::endl;
00680 #endif
00681       stateOnDet(initialTSOS, 
00682                  recHits.front().geographicalId().rawId(),
00683                  initialState);
00684       // Create a new Trajectory Seed    
00685       output[ialgo]->push_back(TrajectorySeed(initialState, recHits, alongMomentum));
00686 #ifdef FAMOS_DEBUG
00687       std::cout << "Trajectory seed created ! " << std::endl;
00688 #endif
00689       break;
00690       // End of the loop over seeding algorithms
00691     }
00692     // End on the loop over simtracks
00693   }
00694 
00695   for ( unsigned ialgo=0; ialgo<seedingAlgo.size(); ++ialgo ) { 
00696     std::auto_ptr<TrajectorySeedCollection> p(output[ialgo]);
00697     e.put(p,seedingAlgo[ialgo]);
00698   }
00699 
00700 }
00701 
00702 // This is a copy of a method in 
00703 // TrackingTools/TrajectoryState/src/TrajectoryStateTransform.cc
00704 // but it does not return a pointer (thus avoiding a memory leak)
00705 // In addition, it's also CPU more efficient, because 
00706 // ts.localError().matrix() is not copied
00707 void 
00708 TrajectorySeedProducer::stateOnDet(const TrajectoryStateOnSurface& ts,
00709                                    unsigned int detid,
00710                                    PTrajectoryStateOnDet& pts) const
00711 {
00712 
00713   const AlgebraicSymMatrix55& m = ts.localError().matrix();
00714   
00715   int dim = 5; 
00716 
00717   float localErrors[15];
00718   int k = 0;
00719   for (int i=0; i<dim; ++i) {
00720     for (int j=0; j<=i; ++j) {
00721       localErrors[k++] = m(i,j);
00722     }
00723   }
00724   int surfaceSide = static_cast<int>(ts.surfaceSide());
00725 
00726   pts = PTrajectoryStateOnDet( ts.localParameters(),
00727                                localErrors, detid,
00728                                surfaceSide);
00729 }
00730 
00731 bool
00732 TrajectorySeedProducer::compatibleWithBeamAxis(GlobalPoint& gpos1, 
00733                                                GlobalPoint& gpos2,
00734                                                double error,
00735                                                bool forward,
00736                                                unsigned algo) const {
00737 
00738   if ( !seedCleaning ) return true;
00739 
00740   // The hits 1 and 2 positions, in HepLorentzVector's
00741   XYZTLorentzVector thePos1(gpos1.x(),gpos1.y(),gpos1.z(),0.);
00742   XYZTLorentzVector thePos2(gpos2.x(),gpos2.y(),gpos2.z(),0.);
00743 #ifdef FAMOS_DEBUG
00744   std::cout << "ThePos1 = " << thePos1 << std::endl;
00745   std::cout << "ThePos2 = " << thePos2 << std::endl;
00746 #endif
00747 
00748 
00749   // Create new particles that pass through the second hit with pT = ptMin 
00750   // and charge = +/-1
00751   
00752   // The momentum direction is by default joining the two hits 
00753   XYZTLorentzVector theMom2 = (thePos2-thePos1);
00754 
00755   // The corresponding RawParticle, with an (irrelevant) electric charge
00756   // (The charge is determined in the next step)
00757   ParticlePropagator myPart(theMom2,thePos2,1.,theFieldMap);
00758 
00762   bool intersect = myPart.propagateToBeamCylinder(thePos1,originRadius[algo]*1.);
00763   if ( !intersect ) return false;
00764 
00765 #ifdef FAMOS_DEBUG
00766   std::cout << "MyPart R = " << myPart.R() << "\t Z = " << myPart.Z() 
00767             << "\t pT = " << myPart.Pt() << std::endl;
00768 #endif
00769 
00770   // Check if the constraints are satisfied
00771   // 1. pT at cylinder with radius originRadius
00772   if ( myPart.Pt() < originpTMin[algo] ) return false;
00773 
00774   // 2. Z compatible with beam spot size
00775   if ( fabs(myPart.Z()-z0) > originHalfLength[algo] ) return false;
00776 
00777   // 3. Z compatible with one of the primary vertices (always the case if no primary vertex)
00778   const reco::VertexCollection* theVertices = vertices[algo];
00779   if (!theVertices) return true;
00780   unsigned nVertices = theVertices->size();
00781   if ( !nVertices || zVertexConstraint[algo] < 0. ) return true;
00782   // Radii of the two hits with respect to the beam spot position
00783   double R1 = std::sqrt ( (thePos1.X()-x0)*(thePos1.X()-x0) 
00784                         + (thePos1.Y()-y0)*(thePos1.Y()-y0) );
00785   double R2 = std::sqrt ( (thePos2.X()-x0)*(thePos2.X()-x0) 
00786                         + (thePos2.Y()-y0)*(thePos2.Y()-y0) );
00787   // Loop on primary vertices
00788   for ( unsigned iv=0; iv<nVertices; ++iv ) { 
00789     // Z position of the primary vertex
00790     double zV = (*theVertices)[iv].z();
00791     // Constraints on the inner hit
00792     double checkRZ1 = forward ?
00793       (thePos1.Z()-zV+zVertexConstraint[algo]) / (thePos2.Z()-zV+zVertexConstraint[algo]) * R2 : 
00794       -zVertexConstraint[algo] + R1/R2*(thePos2.Z()-zV+zVertexConstraint[algo]);
00795     double checkRZ2 = forward ?
00796       (thePos1.Z()-zV-zVertexConstraint[algo])/(thePos2.Z()-zV-zVertexConstraint[algo]) * R2 :
00797       +zVertexConstraint[algo] + R1/R2*(thePos2.Z()-zV-zVertexConstraint[algo]);
00798     double checkRZmin = std::min(checkRZ1,checkRZ2)-3.*error;
00799     double checkRZmax = std::max(checkRZ1,checkRZ2)+3.*error;
00800     // Check if the innerhit is within bounds
00801     bool compat = forward ?
00802       checkRZmin < R1 && R1 < checkRZmax : 
00803       checkRZmin < thePos1.Z()-zV && thePos1.Z()-zV < checkRZmax; 
00804     // If it is, just return ok
00805     if ( compat ) return compat;
00806   }
00807   // Otherwise, return not ok
00808   return false;
00809 
00810 }  
00811