CDF II
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 // String Header First  --
003 //------------------------
004 #include <string>
005 
006 //-----------------------
007 // This Class's Header --
008 //-----------------------
009 #include "TrackingUserMods/SiStripNt.hh"
010 
011 //-------------
012 // C Headers --
013 //-------------
014 #include <assert.h>
015 
016 //---------------
017 // C++ Headers --
018 //---------------
019 #include <iostream>
020 
021 //-------------------------------
022 // Collaborating Class Headers --
023 //-------------------------------
024 #include "AbsEnv/AbsEnv.hh"
025 #include "HepTuple/HepHist1D.h"
026 #include "HepTuple/HepHist2D.h"
027 #include "HepTuple/HepNtuple.h"
028 #include "HepTuple/HepFileManager.h"
029 #include "HepTuple/HepRootFileManager.h"
030 #include "TrackingObjects/Tracks/CdfTrackCutBase.hh"
031 #include <algorithm>
032 #include "ElectronObjects/SimpleExtrapolatedTrack.hh"
033 #include "CalorObjects/Pes2dClusterColl.hh"
034 #include "TrackingObjects/Storable/CdfTrackView.hh"
035 #include "TrackingObjects/Storable/CdfTrackColl.hh"
036 #include "TrackingObjects/Tracks/track_cut.hh"
037 #include "VertexAlg/VertexFitter.hh"
038 #include "GeometryBase/CdfDetector.hh"
039 #include "SiliconGeometry/AbsSiDetectorNode.hh"
040 #include "TrackingSI/PathFinder/SiWaferIntersectionSet.hh"
041 #include "TrackingSI/PathFinder/SmarterSiPathFinder.hh"
042 #include "TrackingObjects/SiData/SiStripSet.hh"
043 #include "TrackingObjects/SiData/SiVrbInfoSet.hh"
044 #include "TrackingObjects/Storable/StorableRun2SiStripSet.hh"
045 
046 #include "GeometryBase/CdfDetector.hh"
047 #include "SiliconGeometry/AbsSiNumerology.hh"
048 #include "SiliconGeometry/HWSiNumerology.hh"
049 #include "SiliconGeometry/CdfHalfLadder.hh"
050 #include "SiliconGeometry/CdfHalfLadderSet.hh"
051 #include "TrackingSI/Integrator/SimpleSiliconIntegrator.hh" 
052 
053 #include "TrackingObjects/SiData/ChipStatus.hh"
054 #include "TrackingObjects/SiData/SiChipKey.hh"
055 
056 #include "SvxDaqObjects/SiPedInfoSet.hh"
057 #include "SvxDaqObjects/MostRecentSiCalibration.hh" 
058 #include "SvxDaqObjects/SiDBInfoSet.hh"
059 #include "SvxDaqObjects/SiStripStatus.hh"
060 
061 //#include "TrackingObjects/SiData/SiStrip.hh"
062 //#include "TrackingObjects/SiData/SiDataRep.hh"
063 
064 #include "SiliconGeometry/SiDigiCode.hh"
065 #include "TrackingObjects/SiData/SiDataShift.hh"
066 
067 #include "SimulationObjects/PropagatedSiParticleColl.hh"
068 
069 #include "TrackingUserHL/Utility/CorrectedHitLocations.hh"
070 
071 #include <map>
072 #include <algorithm>
073 using std::map;
074 
075 //-----------------------------------------------------------------------
076 // Local Macros, Typedefs, Structures, Unions and Forward Declarations --
077 //-----------------------------------------------------------------------
078 
079 static const char rcsid[] = "$Id: SiStripNt.cc,v 1.28 2002/10/18 23:17:38 herndon Exp $";
080 
081 //----------------
082 // Constructors --
083 //----------------
084 SiStripNt::SiStripNt( 
085     const char* const theName, 
086     const char* const theDescription )
087     : HepHistModule( theName, theDescription ),
088       _debug("mydebug",this,true),
089       _missingLayer("missingLayer",this,2),
090       _useCircleFit("useCircleFit",this,false),
091       _inputBanksDescription("InputBanksDescr", this, ""),
092       _ptCut("ptCut", this, 0.2),
093       _chipStat(NULL),_svxDaqInfoSet(NULL),_islDaqInfoSet(NULL),
094       _monteCarlo("monteCarlo", this, false),
095       _intrinsic("intrinsic", this, false),
096       cotbeam("default"),
097       svxbeam("default")
098 {
099   commands()->append(&_debug);
100   commands()->append(&_missingLayer);
101   commands()->append(&_useCircleFit);
102   commands()->append(&_inputBanksDescription);
103   commands()->append(&_ptCut);
104   commands()->append(&_monteCarlo);
105   commands()->append(&_intrinsic);
106   _debug.addDescription("\tDebug mode");
107   _missingLayer.addDescription("\tSpecify the layer to study, offline numbering");
108   _useCircleFit.addDescription("\tFit silicon hits in r-phi (COT constrained curvature)");
109   _inputBanksDescription.addDescription(
110   "                        Specify description of SIXD / ISLD Banks to use.\n\
111                            Note: applies only if InputFromBanks = true." );
112   _ptCut.addDescription("\tTranverse momentum cut for track selection.  Selects tracks above cut.");
113   _monteCarlo.addDescription("\tTrue if using Monte Carlo data.");
114   _intrinsic.addDescription("\tTrue if intrinsic residuals have to be found");
115 }
116 //--------------
117 // Destructor --
118 //--------------
119 SiStripNt::~SiStripNt( )
120 {
121   delete siPathFinder;
122   if (_chipStat) delete _chipStat;
123 
124   delete _svxDaqInfoSet;
125   delete _islDaqInfoSet;
126 }
127 
128 //--------------
129 // Operations --
130 //--------------
131 AppResult SiStripNt::bookHistograms( )
132 {
133   const AbsSiDetectorNode *siDetector = CdfDetector::instance()->getSiDetector();
134   siPathFinder         = new SmarterSiPathFinder(siDetector);  
135   HepFileManager* fm = fileManager();
136   StripNtuple = &fm->ntuple("strips");
137   initNtuple=false;
138 
139   return AppResult::OK;
140 }
141 
142 AppResult SiStripNt::beginRun( AbsEvent* aRun ) {
143   // Get beam spot from beam database
144   //CotBeam cotbeam("default");
145   //cotbeam.loadRun();
146 
147   cotbeam.loadRun();
148   svxbeam.loadRun();
149 
150   // Get the information on the chip settings and integrated ladders
151   if (_chipStat) delete _chipStat;
152   _chipStat = new ChipStatus("ofotl_prd_read",AbsEnv::instance()->runNumber());
153 
154   return AppResult::OK;
155 }
156 
157 AppResult SiStripNt::fillHistograms( AbsEvent* anEvent )
158 {
159   HepFileManager* fm = fileManager();
160   
161   // Prepare the ntuples
162   if (!initNtuple) { // Do it only once
163     initNtuple = true;
164     StripNtuple->columnAt("TRK::event",    &eventnum,     (int) 0);
165     StripNtuple->columnAt("TRK::run",      &runnum,       (int) 0);
166     // tracks
167     StripNtuple->columnAt("TRK::pt",       &pt,           (float) 0);
168     StripNtuple->columnAt("TRK::eta",      &eta,          (float) 0);
169     StripNtuple->columnAt("TRK::curv",     &curvature,    (float) 0);
170     StripNtuple->columnAt("TRK::d0",       &d0,           (float) 0);
171     StripNtuple->columnAt("TRK::phi0",     &phi0,         (float) 0);
172     StripNtuple->columnAt("TRK::z0",       &z0,           (float) 0);
173     StripNtuple->columnAt("TRK::deta",     &deta,         (float) 0);
174     StripNtuple->columnAt("TRK::dcurv",    &dcurvature,   (float) 0);
175     StripNtuple->columnAt("TRK::dd0",      &dd0,          (float) 0);
176     StripNtuple->columnAt("TRK::dphi0",    &dphi0,        (float) 0);
177     StripNtuple->columnAt("TRK::dz0",      &dz0,          (float) 0);
178     StripNtuple->columnAt("TRK::cotnhits", &cotnhits,     (int) 0);
179     StripNtuple->columnAt("TRK::cotnax",   &cotnax,       (int) 0);
180     StripNtuple->columnAt("TRK::cotnst",   &cotnst,       (int) 0);
181     StripNtuple->columnAt("TRK::cotchisq", &cotchisq,     (float) 0);
182     StripNtuple->columnAt("TRK::nphi",     &trnphi,       (int) 0);
183     StripNtuple->columnAt("TRK::nsas",     &trnsas,       (int) 0);
184     StripNtuple->columnAt("TRK::nz",       &trnz,         (int) 0);
185     StripNtuple->columnAt("TRK::beamx",    &beamx,           (float) 0);
186     StripNtuple->columnAt("TRK::beamy",    &beamy,           (float) 0);
187     //StripNtuple->columnAt("TRK::db",       &db,           (float) 0);
188     //StripNtuple->columnAt("TRK::zb",       &zb,           (float) 0);
189     //StripNtuple->columnAt("TRK::mom",      &mom,          (float) 0);
190     StripNtuple->columnAt("TRK::trkq",     &trkq,         (int) 0);
191     StripNtuple->columnAt("TRK::theta",    &tktheta,      (float) 0);
192     
193     StripNtuple->columnAt("HIT::nhit", &nhit, (int) 1);
194     StripNtuple->span("HIT::nhit", 0, MaxHit);
195     StripNtuple->columnArray("HIT::hitq",     1,  (float) 0)   .dimension(MaxHit).index("nhit");
196     StripNtuple->columnArray("HIT::hitlay",   1,  (int) 0)     .dimension(MaxHit).index("nhit");
197     StripNtuple->columnArray("HIT::hitside",  1,  (int) 0)     .dimension(MaxHit).index("nhit");
198     StripNtuple->columnArray("HIT::hitdist",  1,  (float) 0)   .dimension(MaxHit).index("nhit");
199     StripNtuple->columnArray("HIT::hiterr",   1,  (float) 0)   .dimension(MaxHit).index("nhit");
200     StripNtuple->columnArray("HIT::hitnstr",  1,  (int) 0)     .dimension(MaxHit).index("nhit");
201     //StripNtuple->columnArray("HIT::hitglor",  1,  (float) 0)   .dimension(MaxHit).index("nhit");
202     //StripNtuple->columnArray("HIT::hitglop",  1,  (float) 0)   .dimension(MaxHit).index("nhit");
203     //StripNtuple->columnArray("HIT::hitgloz",  1,  (float) 0)   .dimension(MaxHit).index("nhit");
204     //StripNtuple->columnArray("HIT::hitintr",  1,  (float) 0)   .dimension(MaxHit).index("nhit");
205     //StripNtuple->columnArray("HIT::hitintp",  1,  (float) 0)   .dimension(MaxHit).index("nhit");
206     //StripNtuple->columnArray("HIT::hitintz",  1,  (float) 0)   .dimension(MaxHit).index("nhit");
207 
208     StripNtuple->columnAt("INT::stat",     &stat,         (int) 0);
209     StripNtuple->columnAt("INT::bw",       &bandWidth,    (int) 0);
210     StripNtuple->columnAt("INT::thrshold", &threshold,    (int) 0);
211     StripNtuple->columnAt("INT::bestate",  &beState,      (int) 0);
212     StripNtuple->columnAt("INT::sincel1a", &timeSinceL1A, (int) 0);
213     StripNtuple->columnAt("INT::intcogph", &sietap,       (float) 0);
214     StripNtuple->columnAt("INT::intcogz",  &sietaz,       (float) 0);
215     StripNtuple->columnAt("INT::phiint",   &phiint,       (float) 0);
216     StripNtuple->columnAt("INT::philoc",   &philoc,       (float) 0);
217     StripNtuple->columnAt("INT::yloc",     &yloc,         (float) 0);
218     StripNtuple->columnAt("INT::zloc",     &zloc,         (float) 0);
219     //StripNtuple->columnAt("INT::zcent",    &zcent,        (float) 0);
220     StripNtuple->columnAt("INT::xint",     &xint,         (float) 0);  //new
221     StripNtuple->columnAt("INT::yint",     &yint,         (float) 0);  //new
222     StripNtuple->columnAt("INT::zint",     &zint,         (float) 0);
223     StripNtuple->columnAt("INT::xintC",    &xintC,         (float) 0);  //new
224     StripNtuple->columnAt("INT::yintC",    &yintC,         (float) 0);  //new
225     StripNtuple->columnAt("INT::zintC",    &zintC,         (float) 0);
226     StripNtuple->columnAt("INT::entangph", &entangph,     (float) 0);
227     StripNtuple->columnAt("INT::entangz",  &entangz,      (float) 0);
228     StripNtuple->columnAt("INT::rodotphi", &rodotphi,     (int) 0);
229     StripNtuple->columnAt("INT::phiwedge", &phiwedge,     (int) 0);
230     StripNtuple->columnAt("INT::chipp",    &chipidp,        (int) 0);
231     StripNtuple->columnAt("INT::chipz",    &chipidz,        (int) 0);
232     StripNtuple->columnAt("INT::bulkhead", &bulkhead,     (int) 0);
233     //StripNtuple->columnAt("INT::wedgechp", &wedgechp,     (int) 0);
234     //StripNtuple->columnAt("INT::wedgechz", &wedgechz,     (int) 0);
235     StripNtuple->columnAt("INT::layer",    &layer,        (int) 0);
236     //StripNtuple->columnAt("INT::chipidp",  &chipidp,      (int) 0);
237     //StripNtuple->columnAt("INT::chipidz",  &chipidz,      (int) 0);
238     //StripNtuple->columnAt("INT::ladder",   &ladder,       (int) 0);
239     //StripNtuple->columnAt("INT::hlad",     &hlad,         (int) 0);
240     //StripNtuple->columnAt("INT::pcos32",   &pcos32,       (float) 0);
241     StripNtuple->columnAt("INT::overlap",  &overlap,      (int) 0);
242 
243     StripNtuple->columnAt("STP::nstrp", &nstrp, (int)1);
244     StripNtuple->span("STP::nstrp", 0, MaxStrip);
245     StripNtuple->columnArray("STP::strpq",   1, (float)-999) .dimension(MaxStrip).index("nstrp");
246     StripNtuple->columnArray("STP::strpped", 1, (float)-999) .dimension(MaxStrip).index("nstrp");
247     StripNtuple->columnArray("STP::strppedd",1, (float)-999) .dimension(MaxStrip).index("nstrp");
248     StripNtuple->columnArray("STP::strpnoi", 1, (float)-999) .dimension(MaxStrip).index("nstrp");
249     StripNtuple->columnArray("STP::strploc", 1, (float)-999) .dimension(MaxStrip).index("nstrp");
250     StripNtuple->columnArray("STP::strpbad", 1, (int)-999)   .dimension(MaxStrip).index("nstrp");
251     StripNtuple->columnArray("STP::strpnum", 1, (int)-999)   .dimension(MaxStrip).index("nstrp");
252     StripNtuple->columnArray("STP::strprel", 1, (int)-999)   .dimension(MaxStrip).index("nstrp");
253 
254     StripNtuple->columnAt("STZ::nstrz", &nstrz, (int)1);
255     StripNtuple->span("STZ::nstrz", 0, MaxStrip);
256     StripNtuple->columnArray("STZ::strzq",   1, (float)-999) .dimension(MaxStrip).index("nstrz");
257     StripNtuple->columnArray("STZ::strzped", 1, (float)-999) .dimension(MaxStrip).index("nstrz");
258     StripNtuple->columnArray("STZ::strzpedd", 1, (float)-999) .dimension(MaxStrip).index("nstrz");
259     StripNtuple->columnArray("STZ::strznoi", 1, (float)-999) .dimension(MaxStrip).index("nstrz");
260     StripNtuple->columnArray("STZ::strzloc", 1, (float)-999) .dimension(MaxStrip).index("nstrz");
261     StripNtuple->columnArray("STZ::strzbad", 1, (int)-999)   .dimension(MaxStrip).index("nstrz");
262     StripNtuple->columnArray("STZ::strznum", 1, (int)-999)   .dimension(MaxStrip).index("nstrz");
263     StripNtuple->columnArray("STZ::strzrel", 1, (int)-999)   .dimension(MaxStrip).index("nstrz");
264 
265     StripNtuple->columnAt("CLP::nclup", &nclup, (int)1);
266     StripNtuple->span("CLP::nclup", 0, MaxClu);
267     StripNtuple->columnArray("CLP::pclustq",   1, (float) -999).dimension(MaxClu).index("nclup");
268     StripNtuple->columnArray("CLP::pdist",     1, (float) -999).dimension(MaxClu).index("nclup");
269     StripNtuple->columnArray("CLP::pdistC",    1, (float) -999).dimension(MaxClu).index("nclup");
270     StripNtuple->columnArray("CLP::pnstrip",   1, (int)   -999).dimension(MaxClu).index("nclup");
271     StripNtuple->columnArray("CLP::phasbad",   1, (int)   -999).dimension(MaxClu).index("nclup");
272     StripNtuple->columnArray("CLP::pneibad",   1, (int)   -999).dimension(MaxClu).index("nclup");
273     StripNtuple->columnArray("CLP::pavnoise",  1, (float) -999).dimension(MaxClu).index("nclup");
274     StripNtuple->columnArray("CLP::pmaxnois",  1, (float) -999).dimension(MaxClu).index("nclup");
275     StripNtuple->columnArray("CLP::psumnois",  1, (float) -999).dimension(MaxClu).index("nclup");
276     StripNtuple->columnArray("CLP::pneff",     1, (float) -999).dimension(MaxClu).index("nclup");
277     StripNtuple->columnArray("CLP::pmxstrpn",  1, (int)   -999).dimension(MaxClu).index("nclup");
278     StripNtuple->columnArray("CLP::pmxstrpq",  1, (float) -999).dimension(MaxClu).index("nclup");
279     StripNtuple->columnArray("CLP::pmxstrpd",  1, (float) -999).dimension(MaxClu).index("nclup");
280     // new (15/09/02) + 3 lines
281     StripNtuple->columnArray("CLP::pstrfrst",  1, (int) -999).dimension(MaxClu).index("nclup");
282     StripNtuple->columnArray("CLP::pstrlast",  1, (int) -999).dimension(MaxClu).index("nclup");
283     StripNtuple->columnArray("CLP::pstrcntr",  1, (int) -999).dimension(MaxClu).index("nclup");
284 
285     StripNtuple->columnAt("CLZ::ncluz", &ncluz, (int)1);
286     StripNtuple->span("CLZ::ncluz", 0, MaxClu);
287     StripNtuple->columnArray("CLZ::zclustq",   1, (float) -999).dimension(MaxClu).index("ncluz");
288     StripNtuple->columnArray("CLZ::zdist",     1, (float) -999).dimension(MaxClu).index("ncluz");
289     StripNtuple->columnArray("CLZ::zdistC",    1, (float) -999).dimension(MaxClu).index("ncluz");
290     StripNtuple->columnArray("CLZ::znstrip",   1, (int)   -999).dimension(MaxClu).index("ncluz");
291     StripNtuple->columnArray("CLZ::zhasbad",   1, (int)   -999).dimension(MaxClu).index("ncluz");
292     StripNtuple->columnArray("CLZ::zneibad",   1, (int)   -999).dimension(MaxClu).index("ncluz");
293     StripNtuple->columnArray("CLZ::zavnoise",  1, (float) -999).dimension(MaxClu).index("ncluz");
294     StripNtuple->columnArray("CLZ::zmaxnois",  1, (float) -999).dimension(MaxClu).index("ncluz");
295     StripNtuple->columnArray("CLZ::zsumnois",  1, (float) -999).dimension(MaxClu).index("ncluz");
296     StripNtuple->columnArray("CLZ::zneff",     1, (float) -999).dimension(MaxClu).index("ncluz");
297     StripNtuple->columnArray("CLZ::zmxstrpn",  1, (int)   -999).dimension(MaxClu).index("ncluz");
298     StripNtuple->columnArray("CLZ::zmxstrpq",  1, (float) -999).dimension(MaxClu).index("ncluz");
299     StripNtuple->columnArray("CLZ::zmxstrpd",  1, (float) -999).dimension(MaxClu).index("ncluz");
300     // new (15/09/02) + 3 lines
301     StripNtuple->columnArray("CLZ::zstrfrst",  1, (int) -999).dimension(MaxClu).index("ncluz");
302     StripNtuple->columnArray("CLZ::zstrlast",  1, (int) -999).dimension(MaxClu).index("ncluz");
303     StripNtuple->columnArray("CLZ::zstrcntr",  1, (int) -999).dimension(MaxClu).index("ncluz");
304   }
305 
306   // get globals for the event
307   eventnum = AbsEnv::instance()->trigNumber();
308   runnum =   AbsEnv::instance()->runNumber();
309 
310   // set up some silicon things
311   // here we get the iterators needed for the event
312   // clusters:
313   const SiHitSet *myHitList = 0;
314   EventRecord::ConstIterator siHitListIter =  
315     EventRecord::ConstIterator(anEvent,"SiHitSet");
316   ConstHandle<SiHitSet> siHitList;
317   if (!siHitListIter.is_valid()) {
318     return AppResult::ERROR;
319   } 
320   siHitList = siHitListIter;
321   myHitList = siHitList.operator->();    
322   const SiHitSet::collection_type & cluster_set=myHitList->contents();
323 
324   SiVrbInfoSet::DetType detType;
325   if (_svxDaqInfoSet) delete _svxDaqInfoSet;
326   _svxDaqInfoSet = new SiVrbInfoSet();
327   detType = SiVrbInfoSet::SVXII;
328   _svxDaqInfoSet->daqInfoSetUp(detType);
329   if (_islDaqInfoSet) delete _islDaqInfoSet;
330   _islDaqInfoSet = new SiVrbInfoSet();
331   detType = SiVrbInfoSet::ISL;
332   _islDaqInfoSet->daqInfoSetUp(detType);
333   
334   // strips:
335   // Unpack from SVX global banks
336   StorableRun2SiStripSet *rawStripSet = new StorableRun2SiStripSet();
337   rawStripSet->set_description(_inputBanksDescription.value());
338   rawStripSet->setSvxSiVrbInfoSet(_svxDaqInfoSet);
339   rawStripSet->setIslSiVrbInfoSet(_islDaqInfoSet);
340   rawStripSet->loadFromBanks(anEvent);
341 
342   EventRecord::ConstIterator siStripListIter = 
343     EventRecord::ConstIterator(anEvent,"StorableRun2SiStripSet");
344   if (!siStripListIter.is_valid()) {
345     std::cout << "invalid strip list!!" << endl;
346     return AppResult::ERROR;
347   } 
348   const StorableRun2SiStripSet * correctedStripSet;
349   StorableRun2SiStripSet_ch ch( siStripListIter );
350   correctedStripSet = &*ch;
351   
352   //
353   // iterate over the tracks 
354   //
355   CdfTrackView_h hTrackView;
356   CdfTrackView::defTracks( hTrackView );
357   int numfound=0;
358   int numtotal=0;
359   for (CdfTrackView::const_iterator 
360          iTrkLink(hTrackView->contents().begin()); 
361        iTrkLink != hTrackView->contents().end();++iTrkLink) {
362 
363     // is it a silicon track? 
364     const CdfTrack_clnk theTrack = *iTrkLink;
365     numtotal++;
366 
367     // -------------------------------------------
368     // TRACK CUTS:
369     // use tracks with pt>1 and at least 2 si hits
370     // -------------------------------------------
371     if ((theTrack->numSIHits()>1)&&(theTrack->pt()>_ptCut.value())) {
372       numfound++;
373 
374       // -------------------------------------------
375       // FILL TRACK INFORMATION
376       // -------------------------------------------
377       HepVector recPars=theTrack->getHelixFit().getParameters();    
378       HepSymMatrix errPars=theTrack->getHelixFit().getErrorMatrix();
379       
380       pt=theTrack->pt();
381       eta = recPars(1);
382       double theta = atan(1./recPars(1));
383       eta =  (theta>0) ?  log(1./tan(theta/2.)) : -log(1./tan(-theta/2.));
384       curvature=  recPars(2);
385       z0=         recPars(3);
386       d0=         recPars(4);
387       phi0=       recPars(5);
388 
389       // this cause crash:
390       //      deta=       sqrt(errPars(1,1));
391       //      dcurvature= sqrt(errPars(2,2));
392       //      dz0=        sqrt(errPars(3,3));
393       //      dd0=        sqrt(errPars(4,4));
394       //      dphi0=      sqrt(errPars(5,5));
395 
396       deta=       errPars(1,1);
397       dcurvature= errPars(2,2);
398       dz0=        errPars(3,3);
399       dd0=        errPars(4,4);
400       dphi0=      errPars(5,5);
401 
402       // if asked for, replace r phi fit with circle fitter 
403       // with constrained curvature
404       // this overwrites the track helix, and overwrites
405       // curvature, d0, phi0, and their errors
406       // everything else is untouched!
407       // so do not use theTrack variable for trajectory info after here
408       
409       Helix theHelix = theTrack->getTrajectory();
410       if (_useCircleFit.value()) {
411         if (fitCircle(theTrack, theHelix)) {
412           rawStripSet->destroy();
413           return AppResult::OK;
414         }
415         // transfer the two parameters changed:
416         // (don't trust the erros returned by fitCircle)
417         d0 = theHelix.getD0();
418         phi0 = theHelix.getPhi0();
419       }
420 
421       trkq = int(abs(curvature) / curvature + 0.5);
422 
423       // beam position
424       // try svx first, then cot, then just zero it.
425       if (svxbeam.isValid()) {
426         beamx = svxbeam.position(z0).x();
427         beamy = svxbeam.position(z0).y();
428       } else if (cotbeam.isValid()) {
429         beamx = cotbeam.position(z0).x();
430         beamy = cotbeam.position(z0).y();
431       } else {
432         beamx = 0.;
433         beamy = 0.;
434       }
435 
436       // cot info
437       cotnhits=theTrack->numCTHits();
438       cotchisq=theTrack->getHelixFit().getChiSquared()/
439         theTrack->getHelixFit().getNumDegreesOfFreedom();
440       cotnax = theTrack->numCTHitsSt();
441       cotnst = theTrack->numCTHitsAx();
442 
443       // ----------------------------------------------
444       // FILL TRACK HIT INFO
445       // ----------------------------------------------
446       nhit=0;
447       for(CdfTrack::SiHitIterator hitIter = theTrack->beginSIHits();  
448           hitIter != theTrack->endSIHits(); ++hitIter) {
449         // get intersection
450         Trajectory::Location *loc = 
451           theHelix.newIntersectionWith((*hitIter)->getHalfLad()->getPlaneOfWafer());
452         
453         // Corrected hit and intersection
454         CorrectedHitLocations corrHit(*theTrack, **hitIter);
455 
456         if (loc) {
457           HepPoint3D trackIntersection=loc->position();
458           trackIntersection = corrHit.getTrkGlobal();
459           SiDigiCode hitCode = (*hitIter)->getHalfLad()->getDetectorCode();
460           hitside[nhit] = ((*hitIter)->pSideHit()==1) ? 0 : 1;
461           hitlay[nhit]  = (*hitIter)->getSiDigiCode().getLayer();
462           
463           // get perp direction
464           HepVector3D PerpDir;
465           if (hitside[nhit]==0)  {
466             PerpDir  = (*hitIter)->getHalfLad()->getShortDirection();
467           }
468           else {
469             if ((hitlay[nhit]==1)||(hitlay[nhit]==2)||
470                 (hitlay[nhit]==4)) {
471               PerpDir  = (*hitIter)->getHalfLad()->getLongDirection();
472             }
473             else {
474               PerpDir  = (*hitIter)->getHalfLad()->getPerpStereoDirection();
475             }
476           }  
477           hitq[nhit]    =(*hitIter)->getQtotal();
478           hitdist[nhit] =PerpDir.dot((*hitIter)->getGlobalPosition()-trackIntersection);
479           hitdist[nhit] =(corrHit.getHitGlobal() - trackIntersection).mag();
480           hiterr[nhit]  =(*hitIter)->getLocalError();
481           hitnstr[nhit] =(*hitIter)->getNstrip();
482           hitglor[nhit] =(*hitIter)->getGlobalPosition().r();
483           hitglop[nhit] =(*hitIter)->getGlobalPosition().phi();
484           hitgloz[nhit] =(*hitIter)->getGlobalPosition().z();
485           hitintr[nhit] =trackIntersection.r();
486           hitintp[nhit] =trackIntersection.phi();
487           hitintz[nhit] =trackIntersection.z();
488           nhit++;
489         }
490       }
491       
492       // -------------------------------------------------- 
493       // FILL INTERSECTION INFO
494       // -------------------------------------------------- 
495       int num_intersect=0;
496       SiWaferIntersectionSet *intWafers = 
497         siPathFinder->newIntersectionSet(*theTrack,_missingLayer.value()
498                                          ,0.0,0.0);
499       const SiWaferIntersection* intersect = NULL;
500       const CdfHalfLadder* currWafer = NULL;
501       overlap=intWafers->size();
502 
503       // loop over intersections
504       for (SiWaferIntersectionSet::ConstIterator 
505              iwafint=intWafers->begin();
506            iwafint!=intWafers->end(); ++iwafint) {
507         intersect =&(*iwafint);
508         if (intersect->inActiveArea() && 
509             theHelix.newIntersectionWith(intersect->getWafer()->getPlaneOfWafer())
510             ) {
511           num_intersect++;
512 
513           // now we are iterating over the intersections in this layer
514           // in each case, we want: the local coordinate of intersection,
515           // the local coordinate of the nearest strip,
516           // the direction of the intersecting track.
517 
518           // since it hit, get the detector
519           currWafer=intersect->getWafer();
520           SiDigiCode currCode=currWafer->getDetectorCode();
521           
522           // Get the propagated particle
523           const PropagatedSiParticle* sip=0;
524           if( _monteCarlo.value() ){
525             EventRecord::ConstIterator it(anEvent,"PropagatedSiParticleColl");
526             if( it.is_valid() ){
527               ConstHandle<PropagatedSiParticleColl> ph(it);
528               if(! ph->contents().empty() ){
529                 //cout << "*****************" << endl;
530                 for(int i=0; i<ph->contents().size() && !sip; i++){
531                   if( ph->contents()[i].value() == currCode.getKey() ){
532                     sip = &(ph->contents()[i]);
533                     //sip->print();
534                     //cout << "-----------------" << endl;
535                   }
536                 }
537               }
538             }
539             if(! sip){ break; }
540           }
541 
542           // See if this halfladder was integrated in this run or not
543           // and if it was in readall mode
544           map<SiChipKey,Chipstream>::const_iterator ic = 
545             _chipStat->getChipInfo()->find(SiChipKey(currCode,0));
546           if (ic == _chipStat->getChipInfo()->end()) { // NOT integrated
547             stat=-1;
548           }
549           else {
550             stat=0;
551             if (ic->second.readAll() != 0) stat |= 1 << 0;
552             if (ic->second.dpsOn() != 0)   stat |= 1 << 1;
553             bandWidth = ic->second.bandWidth();
554             threshold = ic->second.threshold();
555             if (currCode.getLayer()>0 && currCode.getLayer()<6) {
556               _svxDaqInfoSet->setKey(currCode);
557               beState = _svxDaqInfoSet->beState();
558               timeSinceL1A = _svxDaqInfoSet->vrbInfo()->timeSinceL1A();
559               if (_debug.value()) std::cout << "timeSinceL1A is " << timeSinceL1A << std::endl; 
560               //              _svxDaqInfoSet->reset();
561             }
562             else if (currCode.getLayer()>5){
563               _islDaqInfoSet->setKey(currCode);
564               beState = _islDaqInfoSet->beState();
565               timeSinceL1A = _islDaqInfoSet->vrbInfo()->timeSinceL1A();
566               //              _islDaqInfoSet->reset();
567             } 
568           }         
569 
570           // fill some wafer information about the intersection
571           // getting actual intersection point (use Helix, not Track)
572           // SiWaferInstersection returns the hit strip's location only
573           // also end up with local strip number
574           // repeat everything for phi and z side separately
575           // If "intrinsic" is set, use actual intersection point
576           HepPoint3D globalIntersectionLocation;
577           if( _monteCarlo.value() && _intrinsic.value() ){
578             globalIntersectionLocation =
579               ( sip->entryPoint() + sip->exitPoint() ) / 2;
580           } else {
581             globalIntersectionLocation = theHelix.
582               newIntersectionWith(currWafer->getPlaneOfWafer())->position();
583          } 
584 
585           float localIntersectionPhi=currWafer->
586             getLocalCoord(globalIntersectionLocation,0);
587           float localIntersectionZ=currWafer->
588             getLocalCoord(globalIntersectionLocation,1);
589           // strip numbers intersected
590           strpnint=((int)(currWafer->toStripUnits(localIntersectionPhi,0)+0.5));
591           strznint=((int)(currWafer->toStripUnits(localIntersectionZ,1)+0.5));
592 
593           // check for illegal strip numbers (seems inActiveRegion ain't enough)
594           currWafer->getMultiplexing();
595           if ((strpnint<0)||(strpnint>=currWafer->getPhiStripSpecification()->getNumStrips()))
596             continue;
597           if ((strznint<0)||(strznint>=currWafer->getZStripSpecification()->getNumStrips()*
598                              currWafer->getMultiplexing())) 
599             if (_missingLayer.value()>0) continue;
600 
601           // find where relative to the center of the strip;
602           sietap=1.*strpnint-currWafer->toStripUnits(localIntersectionPhi,0);
603           sietaz=1.*strznint-currWafer->toStripUnits(localIntersectionZ,1);
604 
605           phiint=globalIntersectionLocation.phi();
606           if (phiint<0) phiint+=2*PI_CONST;
607           //new(15/05/02) +2 lines
608           xint=globalIntersectionLocation.x();
609           yint=globalIntersectionLocation.y();
610           zint=globalIntersectionLocation.z();
611 
612           // Get glo int corrected for wafer position and bow
613           CorrectedHitLocations corrInt(*theTrack, *currWafer);
614           xintC = corrInt.getTrkGlobal().x();
615           yintC = corrInt.getTrkGlobal().y();
616           zintC = corrInt.getTrkGlobal().z();
617 
618           zcent=currWafer->getGlobalCenterPosition().z();
619           philoc=globalIntersectionLocation.phi()-currWafer->getGlobalCenterPosition().phi();     
620           if (philoc<-PI_CONST) philoc+=2*PI_CONST;       
621           yloc = localIntersectionPhi;
622           zloc = localIntersectionZ;
623 
624           // rodotphi: tells whether readout is in increasing phi or decreasing phi
625           // for now just use brute force: 
626           // increment strip number by one and see what phi does.
627           // decrement if at edge.
628           float tempphicstrip=globalIntersectionLocation(1);
629           float tempphineistrip;
630           if (strpnint-1. < currWafer->getPhiStripSpecification()->getNumStrips()) {
631             tempphineistrip=currWafer->getGlobalPositionVector(currWafer->
632               toLocalCoord(strpnint+1.0,0,0),unsigned(0)).phi();
633           }    
634           else {
635             tempphineistrip=currWafer->getGlobalPositionVector(currWafer->
636               toLocalCoord(strpnint-1.0,0,0),unsigned(0)).phi();
637           }
638           float deltaphi= tempphicstrip-tempphineistrip;
639           rodotphi= (deltaphi>0.) ? 1 : -1;
640 
641           // getting direction at intersection
642           // since newIntersectionWith doesn't calculate the direction,
643           // use pathlength to rho to force direction calculation
644           HepVector3D intersectionDirection=
645             theHelix.getDirection(theHelix.
646                 getPathLengthAtRhoEquals(globalIntersectionLocation.perp()));
647 
648           // entrance angle stuff.
649           HepVector3D waferNorm=currWafer->getPlaneOfWafer().normal();
650           entangph=waferNorm.phi()-intersectionDirection.phi();
651           entangz=waferNorm.theta()-intersectionDirection.theta();
652           tktheta=intersectionDirection.theta();
653 
654           // intersected detector characteristics
655           int barrel = currCode.getBarrel();
656           int halflad = currCode.getHalfLadder();
657           bulkhead=2*barrel+halflad;
658           layer = currCode.getLayer();
659           phiwedge = currCode.getPhiWedge();
660           
661           // a numbering system for [half]ladders
662           ladder = 60*(bulkhead/2) + 12*(layer-1) + phiwedge + 1;
663           hlad = 60 * bulkhead + 12*(layer-1) + phiwedge + 1;
664 
665           // which chip the intersection occurred in
666           int numPhiChipsInLayer=currWafer->getPhiStripSpecification()->getNumStrips()/128;
667           int numZChipsInLayer=currWafer->getZStripSpecification()->getNumStrips()/128;
668           chipidp=strpnint / 128;
669           chipidz= (_missingLayer.value()>0) ? numPhiChipsInLayer+(strznint/128)%numZChipsInLayer : 0;
670           if (layer==2) chipidp+=4; chipidz+=4;
671           if (layer==3) chipidp+=10;chipidz+=10;
672           if (layer==4) chipidp+=20;chipidz+=20;
673           if (layer==5) chipidp+=30;chipidz+=30;
674           chipidp += hlad*44;
675           chipidz += hlad*44;
676 
677           if (_debug.value()) {
678           // spit out some stuff to make sure it makes sense
679           std::cout << "INTERSECTION QUANTITIES: " << endl;
680           std::cout << "localint, strip:" << strpnint << " " 
681                << " sietap " << sietap << endl;
682           std::cout << "globals phi, z, phiNorm " << phiint << " " << zcent << " " << waferNorm.phi() 
683                << " local phi: " << philoc << endl; 
684           std::cout << "readout dot phi " << rodotphi << " entrance angles, phi, z: " 
685                << entangph << " " << entangz << endl;
686           }     
687 
688           // -------------------------------------------
689           // FILL MISSING CLUSTER INFO
690           // -------------------------------------------
691           // now we'll go through the clusters on this wafer (searched for by digiCode)
692           // retrieve clusters ("hits") on this wafer
693           int numSides=1;
694           if (_missingLayer.value()>0) numSides=2;
695           for (int side=0; side<numSides; side++) {
696             currCode.setPnSide(side);
697             nclu=0;
698 
699             // get perp direction
700             HepVector3D PerpDir;
701             if (side==0)  {
702               PerpDir  = currWafer->getShortDirection();
703             }
704             else {
705               if ((layer==1)||(layer==2)||
706                   (layer==4)) {
707                     PerpDir  = currWafer->getLongDirection();
708                   }
709               else {
710                 PerpDir  = currWafer->getPerpStereoDirection();
711               }
712             }  
713 
714             SiHitSet::const_digi_iterator cdi = cluster_set.find(currCode);
715             if ( cdi != cluster_set.end() ) {
716               // sort cluster by closeness
717               double sortArrayDistance[100];
718               int sortArrayIndex[100];
719               std::vector<double> corrDist;
720               int numhits=(*cdi).second.size();
721               if (numhits>100) {
722                 cout << "too many clusters to sort!!" << endl;
723                 continue;
724               }
725               for ( int i = 0 ; i < numhits ; ++i ) {
726                 if (((*cdi).second[i]->pSideHit()&&(side==0))||
727                     ((*cdi).second[i]->nSideHit()&&(side==1))) {
728                     sortArrayDistance[i]=PerpDir.dot((*cdi).second[i]->getGlobalPosition()-
729                                                      globalIntersectionLocation);
730                     sortArrayIndex[i]=i;
731 
732                     // Get hit locations corrected for wafer position and bow
733                     CorrectedHitLocations corrHit(*theTrack, *(*cdi).second[i]);
734                     corrDist.push_back(corrHit.getTrkGlobal().distance(corrHit.getHitGlobal()));
735                 }
736                 else {
737                   std::cout << "SiStripNT: wrong side hit from SiHitSet_cdi.find" << endl;
738                 }
739               }
740               // now sort them with good ol fashioned bucket sort
741               for (int i=0;i<numhits;i++) {
742                 for (int j=i;j<numhits;j++) {
743                   if (fabs(sortArrayDistance[j])<fabs(sortArrayDistance[i])) {
744                     double tempd=sortArrayDistance[i];
745                     sortArrayDistance[i]=sortArrayDistance[j];
746                     sortArrayDistance[j]=tempd;
747                     int tempi=sortArrayIndex[i];
748                     sortArrayIndex[i]=sortArrayIndex[j];
749                     sortArrayIndex[j]=tempi;
750                   }
751                 }
752               }
753               // Sort the corrected distances too
754               std::sort(corrDist.begin(),corrDist.end());
755 
756               // Bork Bork Bork! 
757               for ( int ii = 0; ii < numhits ; ii++ ) {
758                 SiHit *currHit=(*cdi).second[sortArrayIndex[ii]];
759                 if (nclu>=MaxClu) continue;
760                 
761                 // fill in cluster accessors
762                 clustq[side][nclu]=currHit->getQtotal();
763                 dist[side][nclu]=sortArrayDistance[nclu];
764                 distC[side][nclu]=corrDist[nclu];
765                 nstrip[side][nclu]=currHit->getNstrip();
766 
767 
768                 // fill strip info for clusters
769                 maxstripq[side][nclu]=0.;
770                 avnoise[side][nclu]=0.;
771                 sumnoise[side][nclu]=0.;
772                 neff[side][nclu]=0.;
773 
774                 // access the strips in the cluster...
775                 const SiCluster *currCluster=currHit->getCluster();
776 
777                 //.......... new(15/05/02)....................
778 
779                 clusfrst[side][nclu]= currCluster->firstStrip()->getStripNumber(); 
780                 cluslast[side][nclu]= currCluster->lastStrip()->getStripNumber(); 
781                 //        cluscntr[side][nclu]= currCluster->
782                 //............................................
783 
784                 for (SiCluster::Iterator currStrip=currCluster->begin();
785                      currStrip!=currCluster->end();++currStrip) {
786                   
787                   // debug
788                   if (_debug.value()) {
789                     std::cout << (*currStrip) << endl;
790                   }
791                   
792                   float charge=currStrip->dataShift().scaledToFloat(currStrip->getScaledAdc());
793                   float noise=currStrip->dataShift().scaledToFloat(currStrip->getScaledNoise());
794                   sumnoise[side][nclu]+=noise*noise;
795                   avnoise[side][nclu]+=noise;
796                   
797                   if (charge>maxstripq[side][nclu]) {
798                     maxnoise[side][nclu]=noise;
799                     maxstripq[side][nclu]=charge;
800                     maxstripn[side][nclu]=currStrip->getStripNumber();
801                     float maxstriploc=currWafer->toLocalCoord(maxstripn[side][nclu],side,0);
802                     maxstripdist[side][nclu]=maxstriploc-currHit->getLocalAsMeasured();
803                   }
804                   
805                   neff[side][nclu]+=charge*charge;
806                   
807                   // bad strips/neighbors
808                   hasbad[side][nclu] = (currStrip->badStrip()) ? 1 : 0;
809                   if (ii==0) neibad[side][nclu] = (currStrip->hasBadLowNeighbor()) ? 1 : 0;
810                   if ((!neibad)&&(ii==(nstrip[side][nclu]-1)))
811                     neibad[side][nclu] = (currStrip->hasBadHighNeighbor()) ? 1 : 0;
812                 }
813                 
814                 // process the strip sums
815                 avnoise[side][nclu]/=nstrip[side][nclu];
816                 sumnoise[side][nclu]=sqrt(sumnoise[side][nclu]);
817                 neff[side][nclu]=1./(neff[side][nclu]/(clustq[side][nclu]*clustq[side][nclu]));
818                 
819                 // dump to see how much sense we're making
820                 if (_debug.value()) {
821                   std::cout << "__CLUSTER__\n  side " << side << std::setw(1) << nclu << ":  " 
822                             << "\t clustq      "  << clustq[side][nclu]
823                             << "\t dist        "  << dist[side][nclu] 
824                             << "\t nstrip      "  << nstrip[side][nclu] 
825                             << "\n maxstripq   "  << maxstripq[side][nclu] 
826                             << "\t maxstripn   "  << maxstripn[side][nclu] 
827                             << "\t 1st_stripn  "  << clusfrst[side][nclu] 
828                             << "\t last_stripn "  << cluslast[side][nclu] 
829                             << "\n hasbad      "  << hasbad[side][nclu]
830                             << "\t neibad      "  << neibad[side][nclu] 
831                             << "\n avnoise     "  << avnoise[side][nclu] 
832                             << "\t sumnoise    "  << sumnoise[side][nclu] 
833                             << "\t neff        "  <<  neff[side][nclu] 
834                             << "\t maxstripdist " << maxstripdist[side][nclu] << endl; 
835                   
836                 }
837                 nclu++;
838               }
839               if (_debug.value()) {
840                 cout << "found " << numhits << " clusters in said wafer" << endl;
841               }
842             }
843             if (side==0) nclup=nclu;
844             if (side==1) ncluz=nclu;
845           }
846           
847           // -------------------------------------------
848           // FILL MISSING STRIP INFO
849           // -------------------------------------------
850 
851           // now fill the strip array
852           // get digi_iterator for strips on this wafer (digi_code)
853           for (int side=0;side<numSides;side++) {
854             currCode.setPnSide(side);
855             //      SiStripSet::const_digi_iterator di  = siStripList->contents().find(currCode);
856             SiStripSet::const_digi_iterator di  =  rawStripSet->contents().find(currCode);
857             
858             // resolving mistery of empty blocks...
859             if (_debug.value()) {
860               std::cout << "__STRIP_ARRAY__\n  side " << side 
861                         << "  SIZE: " << rawStripSet->contents().size() << endl;
862             }
863             // loop over strips
864             int nstr=0;
865             //      if (di != siStripList->contents().end() ) {
866             if (di != rawStripSet->contents().end() ) {
867               int numstrips=(*di).second.size();
868 
869               // make a copy of the strip array, using multiplexing
870               // information
871               int mplexFactor = (side==0) ? 1 : currWafer->getMultiplexing();
872               int nChannelsInLadder = (side==0) ? currWafer->getPhiStripSpecification()->getNumStrips() : currWafer->getZStripSpecification()->getNumStrips();
873               std::vector <SiStrip> mplexStripSet;
874               mplexStripSet.clear();
875               mplexStripSet.reserve(mplexFactor*numstrips);
876               for ( int iMult = 0; iMult<mplexFactor; iMult++) {
877                 for ( int i = 0 ; i < numstrips ; ++i ) {
878                   SiStrip stp=(*di).second[i];
879                   int currStripNumber = stp.getStripNumber();
880                   int mplexStripNumber = currStripNumber + nChannelsInLadder*iMult;
881                   stp.setStripNumber(mplexStripNumber);
882                   mplexStripSet.push_back(stp);
883                 }
884               }
885 
886               for ( int i = 0 ; i < mplexStripSet.size() ; ++i ) {
887                 const SiStrip cstp = mplexStripSet[i];
888                 int currStripNumber = cstp.getStripNumber();
889                 int index = (side==0) ? strpnint-currStripNumber+MaxStrip/2
890                   : strznint-currStripNumber+MaxStrip/2;
891                 if (index>=0&&index<MaxStrip) {
892                   str_q[side][nstr]=cstp.dataShift().scaledToFloat(cstp.getScaledAdc());
893                   //              str_noi[side][nstr]=cstp.dataShift().scaledToFloat(cstp.getScaledNoise());
894                   //              str_bad[side][nstr]= (cstp.badStrip()) ? 1 : 0;
895                   str_num[side][nstr]=currStripNumber; 
896                   str_local[side][nstr] = currWafer->toLocalCoord(float(currStripNumber%nChannelsInLadder), 
897                                                        side, int(currStripNumber/nChannelsInLadder));
898                   str_rel[side][nstr]=index;
899 
900                   // pedestals
901                   // yes this is copied from Dr. Maksimovic
902                   const SiDBInfoSet *dbSet=MostRecentSiCalibration::dbInfoSet();
903                   if(dbSet){
904                     SiDBInfoSet::const_digi_iterator dbSet_cdi=  dbSet->contents().find(currCode);
905                     if (dbSet_cdi== dbSet->contents().end()) str_ped[side][nstr]=-1000.;
906                     else {
907                       str_ped[side][nstr]=cstp.dataShift().
908                         scaledToFloat(((*dbSet_cdi).second)[currStripNumber%nChannelsInLadder].pedestal().scaledPedestal(0)); 
909                       str_noi[side][nstr]=cstp.dataShift().
910                         scaledToFloat(((*dbSet_cdi).second)[currStripNumber%nChannelsInLadder].pedestal().scaledNoise(0)); 
911                       str_bad[side][nstr]=((*dbSet_cdi).second)[currStripNumber%nChannelsInLadder].status().isGood() ? 0 : 1;
912                       
913                       // To get the DPS pedestal, subtract the DPS-corrected
914                       // charge from the raw charge.
915                       SiStripSet::ConstIterator corrdi =
916                         correctedStripSet->findStrip(currCode,currStripNumber%nChannelsInLadder);
917                       if (corrdi != correctedStripSet->end()) {
918                         SiStrip corrstrip = *corrdi;
919                         str_pedd[side][nstr] = str_q[side][nstr] -
920                           corrstrip.dataShift().scaledToFloat(corrstrip.getScaledAdc());
921                       }
922                       else {
923                         std::cout << "The raw strip " << currCode << "Stp= " << currStripNumber 
924                                   << " " << cstp << " is not found in the corrected set!" << std::endl;
925                         str_pedd[side][nstr] = -9999.0;
926                       }
927                       
928                       //double noise=(&(*dbSet_cdi).second)->operator[](currStripNumber)
929                       // .pedestal().scaledNoise(0);
930                     }
931                   }
932                   // pedestals (obsolete but easy JB methods)
933                   //str_ped[side][nstr]=-1000;
934                   //if ((pedSet)&&(pedSet->nStrips()>0)) {
935                   //  SiPedInfoSet::ConstIterator pedIt    
936                   //   = pedSet->findStrip(currCode,currStripNumber);
937                   //  if (pedIt != pedSet->end()) {
938                   //    const SiPedInfo pedInfo = *pedIt;
939                   //    str_ped[side][nstr]=pedInfo.mean();
940                   //  }
941                   //double noise=pedInfo.noise(); 
942                   nstr++;
943                 }
944                 if (_debug.value()) {
945                   std::cout << cstp << " curr dist " << index << endl;
946                 }
947               }
948               if (_debug.value()) {
949                 if (numstrips>0) std::cout << "found " << numstrips << " strips in said wafer" << endl;
950               }
951             }
952             if (side==0) nstrp=nstr;
953             if (side==1) nstrz=nstr;
954           }
955           // now fill the ntuple cause we're done
956           StripNtuple->capture("STP::strpq",    &(str_q[0][0]));
957           StripNtuple->capture("STP::strpped",  &(str_ped[0][0]));
958           StripNtuple->capture("STP::strppedd",  &(str_pedd[0][0]));
959           StripNtuple->capture("STP::strpnoi",  &(str_noi[0][0]));
960           StripNtuple->capture("STP::strploc",  &(str_local[0][0]));
961           StripNtuple->capture("STP::strpbad",  &(str_bad[0][0]));
962           StripNtuple->capture("STP::strpnum",  &(str_num[0][0]));
963           StripNtuple->capture("STP::strprel",  &(str_rel[0][0]));
964           StripNtuple->capture("STZ::strzq",    &(str_q[1][0]));
965           StripNtuple->capture("STZ::strzped",  &(str_ped[1][0]));
966           StripNtuple->capture("STZ::strzpedd",  &(str_pedd[1][0]));
967           StripNtuple->capture("STZ::strznoi",  &(str_noi[1][0]));
968           StripNtuple->capture("STZ::strzloc",  &(str_local[1][0]));
969           StripNtuple->capture("STZ::strzbad",  &(str_bad[1][0]));
970           StripNtuple->capture("STZ::strznum",  &(str_num[1][0]));
971           StripNtuple->capture("STZ::strzrel",  &(str_rel[1][0]));
972 
973           StripNtuple->capture("CLP::pclustq",   &(clustq[0][0]));
974           StripNtuple->capture("CLP::pdist",     &(dist[0][0]));
975           StripNtuple->capture("CLP::pdistC",    &(distC[0][0]));
976           StripNtuple->capture("CLP::pnstrip",   &(nstrip[0][0]));
977           StripNtuple->capture("CLP::phasbad",   &(hasbad[0][0]));
978           StripNtuple->capture("CLP::pneibad",   &(neibad[0][0]));
979           StripNtuple->capture("CLP::pavnoise",  &(avnoise[0][0]));
980           StripNtuple->capture("CLP::pmaxnois",  &(maxnoise[0][0]));
981           StripNtuple->capture("CLP::psumnois",  &(sumnoise[0][0]));
982           StripNtuple->capture("CLP::pneff",     &(neff[0][0]));
983           StripNtuple->capture("CLP::pmxstrpn",  &(maxstripn[0][0]));
984           StripNtuple->capture("CLP::pmxstrpq",  &(maxstripq[0][0]));
985           StripNtuple->capture("CLP::pmxstrpd",  &(maxstripdist[0][0]));
986           //new(15/05/02) +3 lines
987           StripNtuple->capture("CLP::pstrfrst",  &(clusfrst[0][0]));
988           StripNtuple->capture("CLP::pstrlast",  &(cluslast[0][0]));
989           // StripNtuple->capture("CLP::pstrcntr",  &(cluscntr[0][0]));
990 
991           StripNtuple->capture("CLZ::zclustq",   &(clustq[1][0]));
992           StripNtuple->capture("CLZ::zdist",     &(dist[1][0]));
993           StripNtuple->capture("CLZ::zdistC",    &(distC[1][0]));
994           StripNtuple->capture("CLZ::znstrip",   &(nstrip[1][0]));
995           StripNtuple->capture("CLZ::zhasbad",   &(hasbad[1][0]));
996           StripNtuple->capture("CLZ::zneibad",   &(neibad[1][0]));
997           StripNtuple->capture("CLZ::zavnoise",  &(avnoise[1][0]));
998           StripNtuple->capture("CLZ::zmaxnois",  &(maxnoise[1][0]));
999           StripNtuple->capture("CLZ::zsumnois",  &(sumnoise[1][0]));
1000           StripNtuple->capture("CLZ::zneff",     &(neff[1][0]));
1001           StripNtuple->capture("CLZ::zmxstrpn",  &(maxstripn[1][0]));
1002           StripNtuple->capture("CLZ::zmxstrpq",  &(maxstripq[1][0]));
1003           StripNtuple->capture("CLZ::zmxstrpd",  &(maxstripdist[1][0]));
1004           //new(15/05/02) +3 lines
1005           StripNtuple->capture("CLZ::zstrfrst",  &(clusfrst[1][0]));
1006           StripNtuple->capture("CLZ::zstrlast",  &(cluslast[1][0]));
1007           // StripNtuple->capture("CLZ::zstrcntr",  &(cluscntr[1][0]));
1008 
1009           StripNtuple->capture("HIT::hitq",      &(hitq[0]));
1010           StripNtuple->capture("HIT::hitlay",    &(hitlay[0]));
1011           StripNtuple->capture("HIT::hitside",   &(hitside[0]));
1012           StripNtuple->capture("HIT::hitdist",   &(hitdist[0]));
1013           StripNtuple->capture("HIT::hiterr",    &(hiterr[0]));
1014           StripNtuple->capture("HIT::hitnstr",   &(hitnstr[0]));
1015           //StripNtuple->capture("HIT::hitglor",   &(hitglor[0]));
1016           //StripNtuple->capture("HIT::hitglop",   &(hitglop[0]));
1017           //StripNtuple->capture("HIT::hitgloz",   &(hitgloz[0]));
1018           //StripNtuple->capture("HIT::hitintr",   &(hitintr[0]));
1019           //StripNtuple->capture("HIT::hitintp",   &(hitintp[0]));
1020           //StripNtuple->capture("HIT::hitintz",   &(hitintz[0]));
1021 
1022           StripNtuple->capture();
1023           StripNtuple->storeCapturedData();
1024         } // in active area
1025       }
1026       if (_debug.value()) {
1027         cout << "Number of intersection: " << num_intersect <<endl;
1028       }
1029     }
1030   }
1031   rawStripSet->destroy();
1032       
1033   if (numtotal&&_debug.value())
1034     std::cout << numfound << " tracks found\n";
1035 
1036   return AppResult::OK;
1037 }
1038 
1039 AppModule*
1040 SiStripNt::clone(const char* cloneName)
1041 {
1042   return new SiStripNt(cloneName,"this module is a clone SiStripNt");
1043 }
1044 
1045 const char * 
1046 SiStripNt::rcsId( ) const
1047 {
1048    return rcsid;
1049 }
1050 
1051 int SiStripNt::fitCircle(const CdfTrack_clnk & theTrack, Helix & theHelix) {
1052   
1053   // get parameters
1054   HepVector recPars=theTrack->getHelixFit().getParameters();    
1055   HepSymMatrix errPars=theTrack->getHelixFit().getErrorMatrix();
1056 
1057   curvature=  recPars(2);
1058   d0=         recPars(4);
1059   phi0=       recPars(5);
1060   dd0=        sqrt(errPars(4,4));
1061   dphi0=      sqrt(errPars(5,5));
1062 
1063   // get into center of circle/radius system 
1064   double R0=abs(0.5/curvature);
1065   HepVector c0(2);
1066   double trkq = abs(curvature) / curvature;
1067   c0[0] = (trkq*d0+R0)*cos(phi0+trkq*3.141593/2.);
1068   c0[1] = (trkq*d0+R0)*sin(phi0+trkq*3.141593/2.);
1069 
1070   // now we have the circle coordinates which form the seed 
1071   // (and curvature constraint),
1072   // get the hits off the track and cache them
1073   double hitx[8],hity[8],hitsx[8],hitsy[8];
1074   int numhits=0;
1075   for (CdfTrack::SiHitIterator hits=theTrack->beginSIHits();
1076        hits!=theTrack->endSIHits();++hits) {
1077     // only use rphi hits of course. 
1078     if ((*hits)->pSideHit()) {
1079       // first get the hit position and error
1080       HepPoint3D globalpos=(*hits)->getGlobalPosition();
1081       HepVector3D  globalerr = (*hits)->pSideHit()   ? 
1082         (*hits)->getHalfLad()->getShortDirection() : 
1083         (*hits)->getHalfLad()->getPerpStereoDirection();
1084       globalerr*=(*hits)->getLocalError();
1085       // cache them
1086       hitx[numhits]=globalpos[0];hity[numhits]=globalpos[1];
1087       hitsx[numhits]=globalerr[0];hitsy[numhits]=globalerr[1];
1088       numhits++;
1089     }
1090   }    
1091 
1092   // iterating over the Newton's method now...
1093   // check for enough hits
1094   if (numhits<2) return(1); 
1095 
1096   HepVector update(2);
1097   HepMatrix Covar(2,2);
1098   int niter=0;
1099   do {
1100     // calculate the first derivative and jacobean 
1101     // d{\chi^2}/dx0 = -2\Sigma(|x_i-x0|-R)\times\frac{x_i-x0}{|x_i-x0|}
1102     // jacobean is a mess
1103     HepVector Deriv(2,0);
1104     HepMatrix Jacob(2,2,0);
1105     HepMatrix I(2,2,1);
1106     double chisq=0;
1107     // iterate over hits
1108     for (int i=0;i<numhits;i++) {
1109       // translate to HepVectors
1110       HepVector x_i(2);x_i[0]=hitx[i];x_i[1]=hity[i];
1111       HepVector s_i(2);s_i[0]=hitsx[i];s_i[1]=hitsy[i];
1112 
1113       // some aux variables for this point:
1114       HepVector Rvec=x_i-c0;
1115       double sigma=(dot(s_i,Rvec))/Rvec.norm();
1116       
1117       // update Deriv,Jacobean. Ignore d(sigma)/d(c0) as a perturbation
1118       // first derivative 
1119       Deriv+= -(2/(sigma*sigma))*(Rvec.norm()-R0)*(Rvec/Rvec.norm());
1120       // Jacobean first term
1121       Jacob += (2/(sigma*sigma))*(Rvec*Rvec.T())*(R0/pow(Rvec.norm(),3));
1122       // Jacobean second term
1123       Jacob -= (2/(sigma*sigma))*(R0/Rvec.norm()-1.)*I;
1124       // chisquare
1125       chisq+=(((Rvec.norm())-R0)*(Rvec.norm()-R0))/(sigma*sigma);
1126     }  
1127     // now we have the sums, we update the circle center
1128     // c0_k+1=c0_k-(Jacob)^-1*Deriv
1129     int err;
1130     Covar = Jacob.inverse(err);
1131     if (err) 
1132       std::cout << "SiStripNt: unable to invert circle fit Jacobean... " << endl;
1133     update = -1.*Covar*Deriv;
1134     c0 += update;
1135     niter++;
1136   } while((update.norm()>0.0001)&&(niter<100));
1137 
1138   if (niter==100) return(1);
1139 
1140   HepVector dc0(2);
1141   dc0[0]=sqrt(0.5*fabs(Covar[0][0]));dc0[1]=sqrt(0.5*fabs(Covar[1][1]));
1142   double varc001=0.5*Covar[0][1];
1143 
1144   // invert to track parameters, overwriting the globals 
1145   phi0=atan2(c0[1],c0[0]);
1146   phi0-=trkq*3.141593/2.;
1147   if (phi0<0.) phi0+=2*3.141592;
1148   d0=trkq*(c0.norm()-R0);
1149 
1150   // and the errors
1151   dd0 =   sqrt(fabs((2*varc001)*((c0[0]*c0[1])/(c0.norm()*c0.norm())) 
1152                + (dc0[0]*dc0[0]*c0[0]*c0[0])/(c0.norm()*c0.norm()) 
1153                + (dc0[1]*dc0[1]*c0[1]*c0[1])/(c0.norm()*c0.norm())));
1154   dphi0 = sqrt(fabs((2*varc001)*((-c0[0]*c0[1])/(pow(c0.norm(),4)))
1155                + (dc0[0]*dc0[0]*c0[0]*c0[0])/(pow(c0.norm(),4))
1156                + (dc0[1]*dc0[1]*c0[1]*c0[1])/(pow(c0.norm(),4))));
1157 
1158   // redefine the helix
1159   theHelix.setD0(d0);
1160   theHelix.setPhi0(phi0);
1161   return(0);
1162 }
1163
[ source navigation ] [ diff markup ] [ identifier search ] [ freetext search ] [ file search ]
This page was automatically generated by the LXR engine.
The LXR team
Valid HTML 4.01!
Send problems or questions to cdfcode@fnal.gov