Project:/TrackingUserMods/src/SiHitDropModule.cc

CDF KITS
[ source navigation ]
[ diff markup ]
[ identifier search ]
[ freetext search ]
[ file search ]
Architecture:	[ i386 ]
Version:	[ 4.10.4 ] [ 4.10.5 ] [ 4.8.4 ] [ 4.8.4l3s ] [ 4.8.5 ] [ 4.9.0 ] [ 4.9.1 ] [ 4.9.1.hpt3 ] [ 4.9.1hpt3 ] [ 4.9.1top1 ] [ 5.0.0 ] [ 5.1.0 ] [ 5.1.0beamonly ] [ 5.1.1 ] [ 5.2.0 ] [ 5.3.0 ] [ 5.3.1 ] [ 5.3.1dsp ] [ 5.3.3 ] [ 5.3.3_nt ] [ 5.3.4 ] [ 6.1.1 ] [ 6.1.1b ] [ 6.1.2 ] [ 6.1.3 ] [ 6.1.4 ] [ 6.1.4int3 ] [ 6.1.4mc ] [ 6.1.4mc_a ] [ 6.1.6 ] [ development ]
001 //---------------------------------------------------------------------------------
002 //SiHitDropModule.hh
003 //SiHitDropModule first copy a track from default trackView, then drop some 
004 //SiHits from the track according to user's reqirements(like drop one layer's 
005 //SiHits) and refit it. This module puts the droped SiHits back to the 
006 //refitted track and append it to the event so that you can pick it up in later 
007 //module. Temperaily this mudole calculated unbiased residual for droped SiHits
008 //and put it in the _residNtuple. You can change your reqirement to drop SiHits 
009 //by a talk-to menu.
010 //Exapmle:
011 //talk SiHitDropModule
012 //drop set t     #choose drop or not,default=true
013 //layer set 3    #specify which layer to drop,default=2
014 //side set 2     #specify which side to drop,0=pSide,1=nSide,2=any ,default=2 
015 //barrel set 9   #specify which barrel to drop,9-do not care which barrel,default=9
016 //exit
017 //-------------------------------------------------------------------------
018 //Kai Yi, et al
019 //yik@pha.jhu.edu
020 //The Johns Hopkin University
021 //
022 //------------------------------------------------------------------------
023 
024 
025 #include "TrackingUserMods/SiHitDropModule.hh"
026 #include "AbsEnv/AbsEnv.hh"
027 #include "GeometryBase/CdfDetector.hh"
028 #include "SiliconGeometry/AbsSiDetectorNode.hh"
029 #include "TrackingSI/ClusterFinding/SiClusterFinder.hh"
030 // #include "TrackingSI/ClusterFinding/SiThreshClusterFinder.hh"
031 #include "TrackingSI/ClusterFinding/SiClusterCorrector.hh"
032 #include "TrackingObjects/SiData/SiClusterSet.hh"
033 #include "TrackingObjects/SiData/SiHitSet.hh"
034 #include "TrackingObjects/Storable/StorableRun2SiStripSet.hh"
035 #include "Edm/ConstHandle.hh"
036 #include "TrackingObjects/Tracks/CdfTrack.hh"
037 #include "TrackingObjects/Storable/CdfTrackColl.hh"
038 #include "TrackingUtils/Integrator/MaterialIntegrator.hh"
039 #include "TrackingSI/TrackFitting/SiTrackFitter.hh"
040 #include "TrackingSI/Integrator/SimpleSiliconIntegrator.hh"
041 #include "HepTuple/HepHBookFileManager.h"
042 #include "HepTuple/HepHBookNtuple.h"
043 // For G3XIntegrator
044 #include "Extrapolator/G3XIntegrator.hh"
045 // For track <--> MC particle matching
046 #include "TrackingUtils/Trybos/TRYGenPrimaryVertexSet.hh"
047 #include "TrackingUtils/GeneratorLevel/GenPrimaryVertex.hh"
048 //#include "TrackingUtils/GeneratorLevel/GenPrimaryVertexSet.hh"
049 #include "TrackingUserHL/Matching/HitContentTPAssociator.hh"
050 
051 //-----------------------------------------------------------------------
052 // Local Macros, Typedefs, Structures, Unions and Forward Declarations --
053 //-----------------------------------------------------------------------
054 
055 //static const char rcsid[] = "$Id: SiHitDropModule.cc,v 1.6 2002/10/18 23:17:38 herndon Exp $";
056 
057 
058 //----------------
059 // Constructors --
060 //----------------
061 SiHitDropModule::SiHitDropModule(const char* const theName, 
062                           const char* const theDescription )
063   : HepHistModule( theName, theDescription ),_drop("drop",this,true),
064     _layer("layer",this,2),_side("side",this,2),_barrel("barrel",this,9),
065     _debug("debug",this,true),_scaleP("scaleP",this,1.0),_matInt("matInt",this,1)
066 {
067   _drop.addDescription("\tThis will let you choose to drop SiHits or not,if not, just do a refit");
068   commands()->append(&_drop);
069   _layer.addDescription("\tThis will let you choose to drop SiHits for which layer if you choose to drop");
070   commands()->append(&_layer);
071   _side.addDescription("\tThis will let you choose to drop SiHits on which side");
072   commands()->append(&_side);
073   _barrel.addDescription("\tThis will let you choose to drop which barrel.9--any barrel");
074   commands()->append(&_barrel);
075   _debug.addDescription("\tThis will let you make debug NTuples, Default: True");
076   commands()->append(&_debug);
077   _scaleP.addDescription("\tThis is the error scaling factor for precursor, Default: 1.0");
078   commands()->append(&_scaleP);
079   _matInt.addDescription("\tMaterial Integrator:  1: Simple,  Else: GEANT, Default: 1");
080   commands()->append(&_matInt);
081   
082 
083 }
084 //--------------
085 // Destructor --
086 //--------------
087 SiHitDropModule::~SiHitDropModule( )
088 {
089 
090 }
091 
092 //--------------
093 // Operations --
094 //--------------
095 
096 AppResult SiHitDropModule:: beginJob(AbsEvent* aJob ) 
097 {
098 
099   HepFileManager* manager = fileManager( );
100   assert( 0 != manager );
101 
102   _ptDiffNtuple=&manager->ntuple("ptDiffNtuple");
103   _residNtuple=&manager->ntuple("residNtuple");
104 
105   //std::cout<<"I am in begin job."<<std::endl;
106   return AppResult::OK;
107 }
108 
109 AppResult SiHitDropModule::beginRun( AbsEvent* aRun )
110 {
111   //std::cout<<"I am in begin run."<<std::endl;
112   return AppResult::OK;
113 }
114 
115 AppResult SiHitDropModule::endRun( AbsEvent* aRun )
116 {
117   //std::cout<<"I am in end run."<<std::endl;
118   return AppResult::OK;
119 }
120 
121 AppResult SiHitDropModule::endJob( AbsEvent* aJob )
122 {
123   //std::cout<<"I am in end job."<<std::endl;
124   return AppResult::OK;
125 }
126 
127 AppResult SiHitDropModule::abortJob( AbsEvent* aJob )
128 {
129   //std::cout<<"I am in abort job."<<std::endl;
130   return AppResult::OK;
131 }
132 
133 
134 AppResult SiHitDropModule::event( AbsEvent* anEvent )
135 {
136   
137   //print the trig number and run number
138   int nEvent = AbsEnv::instance()->trigNumber();
139   int nRun = AbsEnv::instance()->runNumber();
140   
141   if (_debug.value()) {
142     std::cout <<std::endl;
143     std::cout << "SiHitDropModule::trigger # "
144          << std::setw(4) << nRun << "/" << std::setw(5) << nEvent
145          << std::endl;
146   }
147      
148   //define a CdfTrackColl handle 
149   Handle<CdfTrackColl>  refittedTrackSet(new CdfTrackColl() );
150   refittedTrackSet->set_description("AbrigeTracks" );
151 
152   //Loop over all tracks inside defaultTracks
153   CdfTrackView_h view_hndl;
154   CdfTrackView::Error trackStatus=CdfTrackView::defTracks(view_hndl);
155   const CdfTrackView::CollType & tracks=view_hndl->contents();
156   // For track <--> MC association
157   TRYGenPrimaryVertexSet gpvSet(*anEvent);
158   gpvSet.initializeFromObs();
159   if ( ! gpvSet.size()) gpvSet.initializeFromHepevt();
160   TPAssociator *anAssoc = new HitContentTPAssociator(&gpvSet,&*view_hndl);
161 
162   for(CdfTrackView::const_iterator itrack=tracks.begin();
163       itrack!=tracks.end();++itrack)
164     {
165 
166       //Track quality cuts
167       if((*itrack)->numSIHits()==0 || (*itrack)->pt()<1.0 || (*itrack)->numSIHits()<6) continue;
168 
169       // Associated MC Particle  (default: pion)
170       const ParentedParticle *p = NULL;
171       if (anAssoc) p = anAssoc->match(**itrack);
172       int pid = (p) ? p->type() : 211;
173 
174       std::vector<const SiHit *> dropedSIHit;
175       //make a copy for this track
176       CdfTrack * myCopy=new CdfTrack(**itrack);    
177 
178       /*
179       Helix helix=myCopy->getTrajectory(); 
180       std::cout<<endl<<"new track"<<std::endl;
181       std::cout<<"numSiHits="<<myCopy->numSIHits()<<" numCTHits="<<myCopy->numCTHits()<<std::endl;
182       std::cout<<"track id="<<myCopy->id()<<" pt="<<myCopy->pt()<<" 1/2radius="<<myCopy->curvature()<<" moment="<<myCopy->momentum()<<std::endl;
183       if(parent) std::cout<<"parent numSiHits="<<parent->numSIHits()<<" parent numCTHits="<< parent->numCTHits()<<std::endl;      
184       */
185 
186 
187       //find the COT parent track for this track to set precursor
188       CdfTrack_clnk parent=myCopy->parent();
189       if(parent) {   
190         while(parent->parent())
191           {
192             parent=parent->parent();
193             //std::cout<<"parent numSiHits="<<parent->numSIHits()<<" parent numCTHits="<< parent->numCTHits()<<std::endl;
194           }
195         if(parent->numSIHits()!=0 || parent->numCTHits()==0) 
196           std::cout<<"Something strange!"<<std::endl;  
197       }
198 
199       // Initialize Silicon stuff
200       //--------------------------
201       MaterialIntegrator* theMaterialIntegrator;
202       if (_matInt.value()==1) {
203         const AbsSiDetectorNode *siDetector = 
204           CdfDetector::instance()->getSiDetector();
205         theMaterialIntegrator=new SimpleSiliconIntegrator(siDetector);
206       } else {
207         G3XIntegrator* theG3XIntegrator=new G3XIntegrator();
208         theG3XIntegrator->setParticle(pid);
209         theMaterialIntegrator=theG3XIntegrator;
210       }
211       SiTrackFitter siFitter(theMaterialIntegrator);
212       siFitter.setPrecursor(parent.operator->());
213       siFitter.setPrecursorErrorScale(_scaleP.value());
214 
215       //the following ntuple is for comparision between old pt and refitted pt for the same
216       //track(no drop), just for investigation purpose
217       if (siFitter.fit(*myCopy)==TrackFitter::FAILED) {
218         std::cout<<"Track refit Failed!"<<std::endl;
219         continue;
220       }  
221       if (_debug.value()) {
222         _ptDiffNtuple->column("nRun",nRun);
223         _ptDiffNtuple->column("nEvent",nEvent); 
224         _ptDiffNtuple->column("algorithm",(*itrack)->algorithm().value());
225         _ptDiffNtuple->column("pt",(*itrack)->pt());
226         _ptDiffNtuple->column("ptFit",myCopy->pt());
227         _ptDiffNtuple->column("ptDiff",(*itrack)->pt()-myCopy->pt()); 
228         _ptDiffNtuple->column("ptDiffPercent",100.0*((*itrack)->pt()-myCopy->pt())/(*itrack)->pt()); 
229         _ptDiffNtuple->captureThenStore();
230         //end for investigation
231       }
232 
233       //loop to find hits belong to one layer and satisfy other reqiurements and drop it
234       if(_drop.value()) {
235         for (CdfTrack::SiHitIterator h1=(*myCopy).beginSIHits();h1!=(*myCopy).endSIHits();++h1) {
236            if ((*h1)->getHalfLad()->getDetectorCode().getLayer()==_layer.value() 
237                && (_side.value()==2 || _side.value()==(*h1)->nSideHit()) 
238                && (_barrel.value()==9 || _barrel.value()==(*h1)->getHalfLad()->getDetectorCode().getBarrel()) ) 
239              { 
240                /*
241                std::cout<<"layer="<<(*h1)->getHalfLad()->getDetectorCode().getLayer()<<std::endl;
242                std::cout<<"pside="<<(*h1)->pSideHit()<<std::endl;
243                std::cout<<"nside="<<(*h1)->nSideHit()<<std::endl;
244                std::cout<<"barrel="<<(*h1)->getHalfLad()->getDetectorCode().getBarrel()<<std::endl;
245                std::cout<<**h1<<std::endl;
246                */
247                dropedSIHit.push_back(*h1);
248                myCopy->dropHit(*h1);
249              }
250         }
251       }
252 
253 
254       //refit the track after drop some SIHITs
255       if (siFitter.fit(*myCopy)==TrackFitter::FAILED) {
256         std::cout<<"Track refit Failed!"<<std::endl;
257         continue;
258       }
259 
260 
261       /*      
262       std::cout<<"after fit"<<std::endl;          
263       std::cout<<"pt="<<myCopy->pt()<<" 1/2radius="<<myCopy->curvature()<<" moment="<<myCopy->momentum()<<std::endl;     
264       */
265 
266       //set a new Algorithm to this copied track so that you can identify in later module
267       //S.B. Using value 0, i.e. no specified algorithm for the time being
268       myCopy->setAlgorithm( CdfTrack::NoSpecifiedAlg );
269 
270       //the following loop adds droped hits back to the copied track and calculates the residual
271       for(int i=0;i<dropedSIHit.size();i++)
272         {
273           //add the dropped SiHits back to copyed track here
274           myCopy->addHit(dropedSIHit[i]);
275           
276           //std::cout<<i<<"th droped SIHIt: "<<*dropedSIHit[i]<<std::endl;
277 
278           //test to calculate residual
279           double resid;   
280           SiDigiCode key;
281           const SiHit* hit=dropedSIHit[i];
282           key = hit->getHalfLad()->getDetectorCode();
283           // Halflader detector code PnSide has no information
284           if (hit->pSideHit()==1) key.setPnSide(0);
285           if (hit->pSideHit()==0) key.setPnSide(1);
286           // get Perp Direction to strips
287           HepVector3D PerpDir;
288           if (key.getPnSide()==0) 
289             {
290               PerpDir  = hit->getHalfLad()->getShortDirection();
291             }
292           if (key.getPnSide()==1) 
293             {
294               if ((key.getLayer()==1)||(key.getLayer()==2)||
295                   (key.getLayer()==4)) 
296                 {
297                   PerpDir  = hit->getHalfLad()->getLongDirection();
298                 }
299               else 
300                 {
301                   PerpDir  = hit->getHalfLad()->getPerpStereoDirection();
302                 }
303             }  
304           Helix helix=myCopy->getTrajectory();
305           Trajectory::Location *loc = helix.newIntersectionWith(hit->getHalfLad()->getPlaneOfWafer());
306           if ( loc )
307             {
308               HepPoint3D intersectionPoint = loc->position();
309               HepPoint3D hitPoint = hit->getGlobalPosition();
310               resid = intersectionPoint.dot(PerpDir) - hitPoint.dot(PerpDir);
311               delete loc;
312             }
313           else
314             {
315               resid = -999.0;
316               errlog( ELerror, "Helix intersection failed" )
317                 << "@testModule::calcResid"
318                 << "Track helix has no intersection with half-ladder that "
319                 << "has a hit on the track"
320                 << endmsg;
321             }
322           //std::cout<<"resid="<<resid<<std::endl;
323           if (_debug.value()) {
324             _residNtuple->column("nRun",nRun);
325             _residNtuple->column("nEvent",nEvent); 
326             _residNtuple->column("residual",resid*10000.0);
327             _residNtuple->captureThenStore();
328           }
329 
330         } //Dropped SiHits for loop
331       
332       refittedTrackSet->contents().push_back( myCopy );
333 
334     }//track for loop
335   
336   if(refittedTrackSet->contents().size())  anEvent->append(refittedTrackSet);
337   return AppResult::OK;
338 }
[ source navigation ]
[ diff markup ]
[ identifier search ]
[ freetext search ]
[ file search ]
This page was automatically generated by the LXR engine.
The LXR team
Send problems or questions to cdfcode@fnal.gov