Main Page | Modules | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members | Related Pages

UberModule.cxx

Go to the documentation of this file.
00001 
00002 //
00003 // 29.05.2003 -- MAK: Changed NtpSRFitVertex to NtpSRVertex
00004 // uggg!  This doesn't work with R0.18.0
00005 //
00006 
00007 #include "Record/RecArrayAllocator.h"
00008 
00009 #include <set>
00010 #include "TROOT.h"
00011 #include "TObjectTable.h"
00012 #include "MINF_Classes/MINFast.h"
00013 #include "CandData/CandRecord.h"
00014 #include "CandData/CandHeader.h"
00015 #include "CalDetSI/CandCalDetSIHandle.h"
00016 #include "CandDigit/CandDigitListHandle.h"
00017 #include "CandDigit/CandDigitHandle.h"
00018 #include "REROOT_Classes/REROOT_Event.h"
00019 #include "REROOT_Classes/REROOT_FLSDigit.h"
00020 #include "REROOT_Classes/REROOT_StdHep.h"
00021 #include "RerootExodus/RerootExodus.h"
00022 #include "Plex/PlexStripEndId.h"
00023 #include "Plex/PlexSEIdAltLItem.h"
00024 #include "Plex/PlexHandle.h"
00025 #include "Plex/PlexPixelSpotId.h"
00026 #include "UgliGeometry/UgliGeomHandle.h"
00027 #include "Navigation/NavSet.h"
00028 #include "Conventions/CalDigitType.h"
00029 #include "Conventions/CalTimeType.h"
00030 #include "Conventions/Munits.h"
00031 #include "Conventions/ReadoutType.h"
00032 #include "Calibrator/Calibrator.h"
00033 #include "Calibrator/CalTempCalibration.h"
00034 #include "MinosObjectMap/MomNavigator.h"
00035 #include "CandNtupleSR/NtpSRRecord.h"
00036 #include "JobControl/JobCResult.h"
00037 #include "JobControl/JobCModuleRegistry.h" // JOBMODULE registration macro
00038 #include "JobControl/JobCommand.h"
00039 #include "MessageService/MsgService.h"
00040 #include "RecoBase/CandTrackListHandle.h"
00041 #include "RecoBase/CandTrackHandle.h"
00042 #include "RecoBase/CandStripListHandle.h"
00043 #include "RecoBase/CandStripHandle.h"
00044 #include "RecoBase/CandSliceListHandle.h"
00045 #include "RecoBase/CandSliceHandle.h"
00046 #include "RecoBase/CandShowerListHandle.h"
00047 #include "RecoBase/CandShowerHandle.h"
00048 #include "RecoBase/CandEventListHandle.h"
00049 #include "RecoBase/CandEventHandle.h"
00050 #include "RecoBase/LinearFit.h"
00051 #include "RecoBase/PropagationVelocity.h"
00052 #include "RawData/RawRecord.h"
00053 #include "RawData/RawDaqHeader.h"
00054 #include "CandTrackSR/CandTrackSRHandle.h"
00055 #include "CandTrackSR/CandTrackSRListHandle.h"
00056 #include "CandFitTrackSR/CandFitTrackSRHandle.h"
00057 #include "CandNtupleSR/NtpSRTrack.h"
00058 #include "CandNtupleSR/NtpSRFiducial.h"
00059 #include "CandNtupleSR/NtpSRVertex.h"
00060 #include "CandNtupleSR/NtpSREvent.h"
00061 #include "CandNtupleSR/NtpSRShower.h"
00062 #include "RawData/RawChannelId.h"
00063 #include "CalDetSI/Helpers.h"
00064 #include "CalDetDST/UberHit.h"
00065 #include "CalDetDST/UberRecord.h"
00066 #include "CalDetDST/UberRecHeader.h"
00067 #include "CalDetDST/UberMC.h"
00068 #include "CalDetDST/UberModule.h"
00069 
00070 #include "DatabaseInterface/DbiResultPtr.h"
00071 #include "CalDetPID/CalDetBeamMomentum.h"
00072 #include "CalDetPID/CandCalDetPIDHandle.h"
00073 #include "CalDetPID/NtpCalDetPID.h"
00074 
00075 #include "Record/SimSnarlRecord.h"
00076 #include "TParticle.h"
00077 
00078 #include <cmath>
00079 
00080 CVSID("$Id: UberModule.cxx,v 1.27 2007/11/11 09:11:39 rhatcher Exp $");
00081 JOBMODULE(UberModule, "UberModule", "A CalDet ntuple, based on NtupleBase");
00082 
00083 using namespace std;
00084 
00085 
00086 
00087 Int_t UberModule::MakePlaneStripIndex(Int_t plane, Int_t strip)
00088 {
00089    // ought to add range checking
00090    return plane*200 + strip;
00091    
00092 }
00093 // MAK: Feb 15, 2005,
00094 // removed all vestiges of ancient fMC_MIP_Conversion, fMC_PE_Conversion hacks
00095 
00096 
00097 //________________________________________________________________________________
00098 UberModule::UberModule():
00099   fIsMC(kFALSE),
00100   fFirstEvent(kTRUE),
00101   fMCDecision(kFALSE),
00102   fRunNumber(0),
00103   fSubRunNumber(-1),
00104   fStarttime(0),
00105   candrec(0),
00106   fCDLH(0),
00107   fCDSI(0),
00108   fStripUidMap(),
00109   fShowerUidMap(),
00110   fTrackUidMap()
00111 //  fEventUidMap()
00112 {
00113 
00114   MSG("UberModule",Msg::kDebug)<<"In UberModule creator"<<endl;
00115 
00116 }//end UberModule()
00117 //________________________________________________________________________________
00118 
00119 UberModule::~UberModule()
00120 {
00121   MSG("UberModule",Msg::kDebug)<<"In ~UberModule"<<endl;
00122 }//end ~UberModule()
00123 
00124 //________________________________________________________________________________
00125 
00126 
00127 JobCResult UberModule::Get(MomNavigator* mom)
00128 {
00129 
00130   // for debugging
00131   //  RecArrayAllocator& allocator = RecArrayAllocator::Instance();
00132   //  allocator.Print(std::cout);
00133   //  std::cout<<"NUberHit = "<<UberHit::GetNUberHit()<<std::endl;
00134   
00135   MSG("UberModule",Msg::kDebug)<<"In UberModule get"<<endl;
00136 
00137   // Check that mom exists.
00138   assert(mom);  
00139 
00140   //ask mom  for a candidate record
00141   candrec=dynamic_cast<CandRecord*>
00142      (mom->GetFragment("CandRecord","PrimaryCandidateRecord"));
00143   if(candrec==0){
00144     MSG("UberModule",Msg::kError)<<"candrec==0"<<endl;
00145     return JobCResult(JobCResult::kFailed);
00146   }
00147   if(fRunNumber==0){
00148     const CandHeader *ch = candrec->GetCandHeader();
00149     if(ch!=0){
00150       fRunNumber=ch->GetRun();
00151       MSG("UberModule",Msg::kDebug)<<"RUN NUMBER "<<fRunNumber<<endl;
00152     }
00153   }
00154   if(fSubRunNumber<0){
00155      RawRecord *rr = dynamic_cast<RawRecord *>(mom->GetFragment("RawRecord"));
00156      if(rr!=0){
00157         const RawDaqHeader *rdh = dynamic_cast<const RawDaqHeader *>
00158            (rr->GetRawHeader());
00159         if(rdh!=0){
00160            fSubRunNumber=rdh->GetSubRun();
00161            MSG("UberModule",Msg::kInfo)<<"Sub run is "<<fSubRunNumber<<endl;
00162         }
00163      }
00164   }
00165 
00166   fCDLH = dynamic_cast<CandDigitListHandle *>
00167     (candrec->FindCandHandle("CandDigitListHandle","canddigitlist"));
00168   if(fCDLH == 0) {//failure
00169     MSG("UberModule",Msg::kError)<<"fCDLH==0"<<endl;
00170     return JobCResult(JobCResult::kFailed);
00171   }
00172 
00173   //now ask candrecord for a CandCalDetSIHandle
00174   fCDSI = dynamic_cast<CandCalDetSIHandle *>
00175     (candrec->FindCandHandle("CandCalDetSIHandle"));
00176   if(fCDSI == 0) {//failure
00177     MSG("UberModule",Msg::kError)<<"fCDSI==0"<<endl;
00178     return JobCResult(JobCResult::kFailed);
00179   }
00180 
00181   
00182   MSG("UberModule",Msg::kDebug)<<"Returning kPassed"<<endl;
00183   
00184   return JobCResult(JobCResult::kPassed); //success
00185 }//end Get
00186 //________________________________________________________________________________
00187 JobCResult UberModule::Reco(MomNavigator *mom)
00188 {
00189   MSG("UberModule",Msg::kDebug)<<"In UberModule Reco"<<endl;
00190 
00191   // MAK: 12 April 2005
00192   // check if there is a cand record
00193   // if not, abort immediately
00194   // this is also checked by the Get() method
00195   // but if filtering is off so are all bets
00196   if(candrec==0){
00197     MSG("UberModule",Msg::kWarning)<<"No CandRecord Found, skipping this routine."<<endl;
00198     return JobCResult(JobCResult::kFailed);
00199   }
00200   
00201   // MAK: 12 April 2005 : originally got the canddigithandleitr hi here
00202   /*
00203   //get an iterator over the digits
00204   CandDigitHandleItr hi(fCDLH->GetDaughterIterator());
00205   if(!hi.IsValid()){
00206     MSG("UberModule",Msg::kError)<<"DaugherIter not valid"<<endl;
00207     return JobCResult(JobCResult::kFailed);
00208   }
00209   */
00210 
00211   //Get a validity context and a plex handle
00212   //
00213   // MAK: 12 April, 2005 : Take VldContext from header
00214   //                       useful if event has no digits
00215   //  const VldContext *vc = (*hi)->GetVldContext();
00216 
00217   VldContext vc_obj = candrec->GetCandHeader()->GetVldContext();
00218   const VldContext* vc = &vc_obj; // for compatibility with much code below
00219   PlexHandle ph(*vc, kTRUE);
00220 
00221   if(!fMCDecision){
00222     if(vc->GetSimFlag()==SimFlag::kReroot||vc->GetSimFlag()==SimFlag::kMC){
00223       fIsMC = kTRUE;
00224       MSG("UberModule",Msg::kInfo)<<"This is an MC file with simflag: "<<SimFlag::AsString(vc->GetSimFlag())<<endl;
00225     }
00226     fMCDecision = kTRUE;
00227   }
00228 
00229   //ReInit the mip calibrator
00230   Calibrator::Instance().Reset(*vc);
00231 
00232   // Use validity context to create header for ntuple record
00233   UberRecHeader newHdr(*vc);
00234   newHdr.SetRunNo(fRunNumber);
00235   newHdr.SetSubRunNo(fSubRunNumber);
00236   
00237   if(fFirstEvent){
00238     fStarttime = vc->GetTimeStamp().GetSec();
00239 
00240     fFirstEvent=kFALSE;
00241   }
00242   newHdr.SetStartTime(fStarttime);
00243 
00244   // Read temperature from the db. NOTE: this comes from RS probe, not DCS!!!!!
00245   float temp = -999.9; // obviously fake value
00246 //  cout<<"about to get temperature"<<endl;
00247 
00248   DbiResultPtr<CalTempCalibration> rptc(*vc);
00249   if(rptc.GetNumRows()>0){
00250      const CalTempCalibration *tc = rptc.GetRow(0);
00251      temp=tc->GetTemp();
00252   }
00253 
00254   if(fIsMC){
00255        // set a bogus temperature
00256     newHdr.SetTemperature(temp);
00257   }
00258   else{
00259     newHdr.SetTemperature(Calibrator::Instance().GetTemperature());
00260   }
00261 
00262   // Read momentum from the db!
00263   float p = 0.1;// obviously fake value
00264   if(fIsMC){
00265 
00266        // try to figure out the momentum from mc truth
00267        // grab the simrecord
00268        const SimSnarlRecord* simrec=dynamic_cast<const SimSnarlRecord*>
00269             (mom->GetFragment("SimSnarlRecord"));
00270        if(simrec){
00271           MSG("UberModule",Msg::kDebug)<<"Got a SimSnarlRecord"<<endl;
00272             // grab the stdhep block
00273             const TClonesArray* stdhep = dynamic_cast<const TClonesArray*> 
00274                  (simrec->FindComponent("TClonesArray", "StdHep"));
00275             if(stdhep){
00276                MSG("UberModule",Msg::kDebug)<<"Got a StdHep"<<endl;
00277                  // primary should be first in the list
00278                  const TParticle* part=
00279                       dynamic_cast<const TParticle*>(stdhep->At(0));
00280                  if(part) {
00281                     MSG("UberModule",Msg::kDebug)<<"Got a Particle: "
00282                                                  <<"  "<<part->GetName()
00283                                                  <<", p= "<<part->P()
00284                                                  <<", id="<<part->GetPdgCode()
00285                                                  <<endl;
00286                     p = part->P();
00287                  }
00288             }
00289        }
00290        else {
00291           MSG("UberModule",Msg::kWarning)
00292              <<"Couldn't get a SimSnarlRecord for this event!"<<endl;
00293        }
00294   }
00295   else{
00296        DbiResultPtr<CalDetBeamMomentum> cdbm(*vc);
00297        if(cdbm.GetNumRows()==1) p = cdbm.GetRow(0)->GetMomentum();
00298   }
00299   newHdr.SetBeamMomentum(p);
00300 
00301   //create a new uberrecord
00302   UberRecord *snarldata = new UberRecord(newHdr);
00303  
00304 //  MSG("UberModule",Msg::kDebug)<<"Resetting Event"<<endl;
00305 //  snarldata->ResetEvent();
00306   fStripUidMap.clear();
00307   fShowerUidMap.clear();
00308   fTrackUidMap.clear();
00309 //  fEventUidMap.clear();
00310 
00311   //start filling the snarldata
00312 
00313   //get the time since beginning of timeframe from from candigitlisthandle
00314   snarldata->triggertime = fCDLH->GetAbsTime();
00315 
00316   //get pid info from the CandCalDetSIHandle
00317   snarldata->snarlno = fCDSI->GetSnarl();
00318   snarldata->triggerword = fCDSI->GetTrigSource();
00319   snarldata->ceradc[0]=fCDSI->GetKovADC2();
00320   snarldata->ceradc[1]=fCDSI->GetKovADC1();
00321   snarldata->ceradc[2]=fCDSI->GetKovADC3();
00322   snarldata->torbits=fCDSI->GetTriggerORBits();
00323   snarldata->torok=fCDSI->GetTriggerOROK();
00324   // MAK: 12 April 2005
00325   // moved here from below
00326   this->FillNtpCalDetPID(snarldata);
00327 
00328 
00329   MSG("UberModule",Msg::kDebug)<<"snarl no "<<snarldata->snarlno<<endl;
00330 
00331   //get time relative to trigger time of cerenkov hits
00332   //note:  CalDetConstants::FARTIMECONVERT()--magic number to convert ticks to 
00333   //nanoseconds for the FAR DETECTOR electronics
00334   //we will need to fix this if the cerenkov is ever read out
00335   //using Near detector electronics (ie 2003)
00336   if(fCDSI->GetKovTimeStamp2()!=0){
00337     snarldata->certime[0]=fCDSI->GetKovTimeStamp2()
00338       *CalDetConstants::FARTIMECONVERT-snarldata->triggertime/Munits::nanosecond;
00339     MSG("UberModule",Msg::kDebug)<<"Kov 2 timestamp "<<fCDSI->GetKovTimeStamp2()
00340                               <<" time convert "<<CalDetConstants::FARTIMECONVERT
00341                               <<" trigger time "<<snarldata->triggertime
00342                               <<" in ns "
00343                               <<snarldata->triggertime/Munits::nanosecond
00344                               <<endl;
00345   }
00346   if(fCDSI->GetKovTimeStamp1()!=0){
00347     snarldata->certime[1]=fCDSI->GetKovTimeStamp1()
00348       *CalDetConstants::FARTIMECONVERT-snarldata->triggertime/Munits::nanosecond;
00349   }
00350   if(fCDSI->GetKovTimeStamp3()!=0){
00351     snarldata->certime[2]=fCDSI->GetKovTimeStamp3()
00352       *CalDetConstants::FARTIMECONVERT-snarldata->triggertime/Munits::nanosecond;
00353   }
00354   snarldata->toftdc[0]=fCDSI->GetTofTDC0();
00355   snarldata->toftdc[1]=fCDSI->GetTofTDC1();
00356   snarldata->toftdc[2]=fCDSI->GetTofTDC2();
00357   snarldata->tofadc[0]=fCDSI->GetTofADC0();
00358   snarldata->tofadc[1]=fCDSI->GetTofADC1();
00359   snarldata->tofadc[2]=fCDSI->GetTofADC2();
00360   if(fCDSI->GetTofTimeStamp()!=0){
00361     snarldata->toftime=fCDSI->GetTofTimeStamp()
00362       *CalDetConstants::FARTIMECONVERT-snarldata->triggertime/Munits::nanosecond;
00363         MSG("UberModule",Msg::kDebug)<<"tof timestamp "<<fCDSI->GetTofTimeStamp()
00364                               <<" time convert "<<CalDetConstants::FARTIMECONVERT
00365                               <<" trigger time "<<snarldata->triggertime
00366                               <<" in ns "
00367                               <<snarldata->triggertime/Munits::nanosecond
00368                                   <<" tof time "<<snarldata->toftime<<endl;
00369 
00370   }
00371   if(fCDSI->GetTofADCTimeStamp0()!=0){
00372     snarldata->tofhittime[0]=fCDSI->GetTofADCTimeStamp0()
00373       *CalDetConstants::FARTIMECONVERT-snarldata->triggertime/Munits::nanosecond;
00374   }
00375   if(fCDSI->GetTofADCTimeStamp1()!=0){
00376     snarldata->tofhittime[1]=fCDSI->GetTofADCTimeStamp1()
00377       *CalDetConstants::FARTIMECONVERT-snarldata->triggertime/Munits::nanosecond;
00378   }
00379   if(fCDSI->GetTofADCTimeStamp2()!=0){
00380     snarldata->tofhittime[2]=fCDSI->GetTofADCTimeStamp2()
00381       *CalDetConstants::FARTIMECONVERT-snarldata->triggertime/Munits::nanosecond;
00382   }
00383 
00384   //add the number of seconds since the beginning of the run to the
00385   //trigger time to make trigger time absolute since beginning of run
00386   snarldata->triggertime+=vc->GetTimeStamp().GetSec()-fStarttime;
00387   MSG("UberModule",Msg::kDebug)<<"Trigger time "<<snarldata->triggertime<<endl
00388                             <<"vld context time "
00389                             <<vc->GetTimeStamp().GetSec()<<endl;
00390 /*
00391   PlexHandle ph(*vc, kTRUE);
00392 
00393   if(!fMCDecision){
00394     if(vc->GetSimFlag()==SimFlag::kReroot||vc->GetSimFlag()==SimFlag::kMC){
00395       fIsMC = kTRUE;
00396     }
00397     fMCDecision = kTRUE;
00398   }
00399 */
00400 
00401 
00402   Bool_t gotsigcorrconv = kFALSE;
00403 
00404 
00405   //check dead chip list in the CandCalDetSIHandle to find mindeadplane
00406   //and and number of dead planes
00407   set<UShort_t> deadplanemap;
00408   vector<RawChannelId>::const_iterator dcit(fCDSI->GetDeadChips().begin());
00409   MSG("UberModule", Msg::kDebug)<<"Starting to get dead chip vector"<<endl
00410                                <<"Size of Dead Chip vector "
00411                                <<fCDSI->GetDeadChips().size()<<endl;
00412   while(dcit!=fCDSI->GetDeadChips().end()){
00413     //make a rawchannel id out of the dead chip (use va channel 3)
00414 //changed 8-27-03 to make dead chip checking work with QIE chips
00415     RawChannelId newid;
00416     newid.SetDetector((*dcit).GetDetector());
00417     newid.SetElecType((*dcit).GetElecType()); 
00418     newid.SetCrate((*dcit).GetCrate());
00419     if((*dcit).GetElecType()==ElecType::kVA){
00420        newid.SetVarcId((*dcit).GetVarcId());
00421        newid.SetVmm((*dcit).GetVmm());
00422        newid.SetVaAdcSel((*dcit).GetVaAdcSel());
00423        newid.SetVaChip((*dcit).GetVaChip());
00424        newid.SetVaChannel(3);
00425     }
00426     if((*dcit).GetElecType()==ElecType::kQIE){
00427        newid.SetGeographicAddress((*dcit).GetGeographicAddress());
00428        newid.SetMasterChannel((*dcit).GetMasterChannel());
00429        newid.SetMinderChannel((*dcit).GetMinderChannel());
00430     }
00431     //setmodebits no longer exists (12/08/04), doe each one separately
00432 //    newid.SetModeBits(false,false,false);
00433     newid.SetPedMode(false);
00434     newid.SetSparsMode(false);
00435     newid.SetCommonMode(false);
00436        
00437     if(ph.GetSEIdAltL(newid).GetPlane()>=0){
00438        deadplanemap.insert(ph.GetSEIdAltL(newid).GetPlane());
00439     }
00440     if(newid.GetElecType()==ElecType::kVA){
00441        if(fCDSI->GetCerenkovChannel3().GetElecType()==ElecType::kVA){
00442           if(newid.IsSameVAChip(fCDSI->GetCerenkovChannel3())){
00443              deadplanemap.insert(65);
00444           }
00445        }
00446        if(fCDSI->GetCerenkovChannel1().GetElecType()==ElecType::kVA){
00447           if(newid.IsSameVAChip(fCDSI->GetCerenkovChannel1())){
00448              deadplanemap.insert(66);
00449           }
00450        }
00451        if(fCDSI->GetCerenkovChannel2().GetElecType()==ElecType::kVA){
00452           if(newid.IsSameVAChip(fCDSI->GetCerenkovChannel2())){
00453              deadplanemap.insert(67);
00454           }
00455        }
00456     }
00457     else if(newid.GetElecType()==ElecType::kQIE){
00458        if(fCDSI->GetCerenkovChannel3().GetElecType()==ElecType::kQIE){
00459           if(newid.IsSameChannel(fCDSI->GetCerenkovChannel3())){
00460              deadplanemap.insert(65);
00461           }
00462        }
00463        if(fCDSI->GetCerenkovChannel1().GetElecType()==ElecType::kQIE){
00464           if(newid.IsSameChannel(fCDSI->GetCerenkovChannel1())){
00465              deadplanemap.insert(66);
00466           }
00467        }
00468        if(fCDSI->GetCerenkovChannel2().GetElecType()==ElecType::kQIE){
00469           if(newid.IsSameChannel(fCDSI->GetCerenkovChannel2())){
00470              deadplanemap.insert(67);
00471           }
00472        }
00473     }
00474 
00475     //make another possible rawchannel id out of the dead chip(use va channel 11)
00476     //which could possibly be on the next plane
00477     //modified 8-17-03 to work with qie elec., 
00478     //need to look at 4 different masterchans to make sure we get all the planes
00479     RawChannelId newid2;
00480     newid2.SetDetector((*dcit).GetDetector());
00481     newid2.SetElecType((*dcit).GetElecType()); 
00482     newid2.SetCrate((*dcit).GetCrate());
00483     if((*dcit).GetElecType()==ElecType::kVA){
00484        newid2.SetVarcId((*dcit).GetVarcId());
00485        newid2.SetVmm((*dcit).GetVmm());
00486        newid2.SetVaAdcSel((*dcit).GetVaAdcSel());
00487        newid2.SetVaChip((*dcit).GetVaChip());
00488        newid2.SetVaChannel(11);
00489     }
00490     if((*dcit).GetElecType()==ElecType::kQIE){
00491        newid2.SetGeographicAddress((*dcit).GetGeographicAddress());
00492        newid2.SetMasterChannel((*dcit).GetMasterChannel()+1);
00493        newid2.SetMinderChannel((*dcit).GetMinderChannel());
00494     }
00495     //setmodebits has disappeared 12-08-04
00496 //    newid2.SetModeBits(false,false,false);
00497     newid2.SetPedMode(false);
00498     newid2.SetSparsMode(false);
00499     newid2.SetCommonMode(false);
00500        
00501 
00502     
00503     if(ph.GetSEIdAltL(newid2).GetPlane()>=0){
00504       deadplanemap.insert(ph.GetSEIdAltL(newid2).GetPlane());
00505     }
00506     if((*dcit).GetElecType()==ElecType::kQIE){
00507        RawChannelId newid3;
00508        newid3.SetDetector((*dcit).GetDetector());
00509        newid3.SetElecType((*dcit).GetElecType()); 
00510        newid3.SetCrate((*dcit).GetCrate());
00511        newid3.SetElecType((*dcit).GetElecType()); 
00512        newid3.SetGeographicAddress((*dcit).GetGeographicAddress());
00513        newid3.SetMasterChannel((*dcit).GetMasterChannel()+2);
00514        newid3.SetMinderChannel((*dcit).GetMinderChannel());
00515 //       newid3.SetModeBits(false,false,false);
00516        newid3.SetPedMode(false);
00517        newid3.SetSparsMode(false);
00518        newid3.SetCommonMode(false);
00519        
00520     
00521        if(ph.GetSEIdAltL(newid3).GetPlane()>=0){
00522           deadplanemap.insert(ph.GetSEIdAltL(newid3).GetPlane());
00523        }
00524        RawChannelId newid4;
00525        newid4.SetDetector((*dcit).GetDetector());
00526        newid4.SetElecType((*dcit).GetElecType()); 
00527        newid4.SetCrate((*dcit).GetCrate());
00528        newid4.SetElecType((*dcit).GetElecType()); 
00529        newid4.SetGeographicAddress((*dcit).GetGeographicAddress());
00530        newid4.SetMasterChannel((*dcit).GetMasterChannel()+3);
00531        newid4.SetMinderChannel((*dcit).GetMinderChannel());
00532 //       newid4.SetModeBits(false,false,false);
00533        newid4.SetPedMode(false);
00534        newid4.SetSparsMode(false);
00535        newid4.SetCommonMode(false);
00536        
00537   
00538        if(ph.GetSEIdAltL(newid4).GetPlane()>=0){
00539           deadplanemap.insert(ph.GetSEIdAltL(newid4).GetPlane());
00540        }
00541     }
00542 
00543 
00544     dcit++;
00545   }
00546   if(deadplanemap.size()!=0){
00547     snarldata->mindeadplaneno = *(deadplanemap.begin());
00548     snarldata->ndeadplanes = deadplanemap.size();
00549   }
00550   MSG("UberModule",Msg::kDebug)<<"First dead plane "<<snarldata->mindeadplaneno
00551                             <<" No. dead planes "<<snarldata->ndeadplanes<<endl;
00552   
00553   //get an iterator over the digits
00554   CandDigitHandleItr hi(fCDLH->GetDaughterIterator());
00555   if(!hi.IsValid()){
00556     // MAK: 12 April 2005
00557     // we can have the case where there are no digits in the event
00558     // it's not common but can happen, particularly for low energy protons
00559     // we'd like to create an uberrecord to record such cases
00560     // ... similar to the case of an ND spill w/o neutrino interactions
00561     MSG("UberModule",Msg::kInfo)<<"DaugherIter not valid: probably no hits in event"<<endl;
00562     MSG("UberModule",Msg::kDebug)<<"giving fragment to mom"<<endl;  
00563     mom->AdoptFragment(snarldata);
00564 
00565     return JobCResult(JobCResult::kAOK);
00566   }
00567 
00568   //sort hits in event by plane
00569   CandDigitHandleKeyFunc* kf = hi.CreateKeyFunc();
00570   kf->SetFun(CalHelpers::KeyFromPlane);
00571   hi.GetSet()->AdoptSortKeyFunc(kf,kTRUE,kFALSE);
00572   //further sort hits by strip
00573   kf=hi.CreateKeyFunc();
00574   kf->SetFun(CalHelpers::KeyFromStrip);
00575   hi.GetSet()->AdoptSortKeyFunc(kf,kFALSE,kFALSE);
00576   //and finally, sort by strip end
00577   kf=hi.CreateKeyFunc();
00578   kf->SetFun(CalHelpers::KeyFromEnd);
00579   hi.GetSet()->AdoptSortKeyFunc(kf,kFALSE,kTRUE);
00580   kf = 0;
00581   MSG("UberModule",Msg::kDebug)<<"Done sorting"<<endl;
00582   
00583   //declare variables to find maxes
00584   Float_t maxplanemip = 0.;
00585   Float_t maxp0stripmip=0.;
00586   Float_t maxp1stripmip=0.;
00587   Float_t maxp0hitmip=0.;
00588   Float_t maxp1hitmip=0.;
00589   Float_t totmipeven=0.;
00590   Float_t totmipodd=0.;
00591 
00592 /*
00593   map<PlexStripEndId, REROOT_FLSDigit> FLSDigitMap;
00594   int nmchitstrips=0;
00595   if(fIsMC){
00596     const TClonesArray *hep = RerootExodus::GetStdHepList();
00597     mc->SetStdHep(hep);
00598     MSG("UberModule",Msg::kDebug)<<"size of hep is "<<hep->GetEntries()<<endl;
00599     REROOT_Event *revt = gMINFast->GetREROOTEvent();
00600     mc->nhep = revt->n_stdheps();
00601     MSG("UberModule",Msg::kDebug)<<"and mc nhep is  "<<mc->nhep<<endl;
00602     REROOT_StdHep *stdhep = static_cast<REROOT_StdHep *>(hep->At(0));
00603     if(stdhep!=0){
00604       mc->mctype = stdhep->ID();
00605       mc->mcpx = stdhep->Px();
00606       mc->mcpy = stdhep->Py();
00607       mc->mcpz = stdhep->Pz();
00608       mc->mcenergy = stdhep->E();
00609       mc->mcvx = stdhep->Xmm();
00610       mc->mcvy = stdhep->Ymm();
00611       mc->mcvz = stdhep->Zmm();
00612     }
00613     const TClonesArray *flsd = RerootExodus::GetFLSDigitList();
00614     int i=0;
00615     TObject *obj;
00616     while((obj=flsd->At(i))){
00617       REROOT_FLSDigit *mcdig = static_cast<REROOT_FLSDigit*>(obj);
00618       PlexStripEndId mcseid = RerootExodus::PECAB2SEId(mcdig->IPln(),
00619                                                        mcdig->IExtr(),
00620                                                        mcdig->ICell(),0);
00621       FLSDigitMap[mcseid]=*mcdig;
00622       i++;
00623     }
00624   }
00625 */
00626     
00627 
00628   
00629   //begin looping over planes
00630   MSG("UberModule",Msg::kDebug)<<"begining to loop over hits"<<endl;;
00631   for(hi.Reset();hi.IsValid();hi.NextKey()){ 
00632     UShort_t plane = (*hi)->GetPlexSEIdAltL().GetPlane();
00633     MSG("UberModule",Msg::kDebug)<<"PLANE "<<plane<<endl;
00634     Bool_t iscosmic = kFALSE;
00635     if(plane>=CalDetConstants::PLANECONST){ //if not a detector plane, then it's a cosmic (bad in 2003)
00636       iscosmic=kTRUE;
00637     }
00638     if(!iscosmic){
00639       snarldata->nhitplanes++; //increment nhitplanes
00640     }
00641     Float_t planemip=0.;
00642 
00643     //begin looping over strips
00644     for(CandDigitHandleItr siter(hi, kTRUE);
00645         siter.IsValid();siter.NextKey()){ 
00646       UShort_t strip = (*siter)->GetPlexSEIdAltL().GetBestSEId().GetStrip();
00647       MSG("UberModule",Msg::kDebug)<<"STRIP "<<strip<<endl;
00648       if(!iscosmic){
00649         MSG("UberModule",Msg::kDebug)<<"adding one to nhitstrips"<<endl;
00650         snarldata->nhitstrips++;  //increment nhitstrips
00651         MSG("UberModule",Msg::kDebug)<<"adding next hit to tclonesarray"<<endl;
00652         //make a hit in UberEvent's TClonesArray
00653         int sindx = snarldata->AddNextHit(plane,strip);
00654         MSG("UberModule",Msg::kDebug)<<"added next hit to tclonesarray"<<endl;
00655         //int uid = (*siter)->GetUidInt();
00656         Int_t uid = MakePlaneStripIndex(plane, strip);
00657         fStripUidMap.insert(std::make_pair(uid,sindx));
00658 
00659 /*
00660         if(fIsMC){
00661           mc->AddNextHit(plane,strip);
00662         }
00663 */
00664       }
00665 /*
00666       if(iscosmic){
00667         snarldata->AddNextCosmicHit(plane,strip);
00668       }
00669 */
00670       Float_t p0stripmip=0.;
00671       Float_t p1stripmip=0.;
00672 
00673       //begin looping over ends
00674       for(CandDigitHandleItr eiter(siter,kTRUE);
00675           eiter.IsValid();eiter.Next()){ //loop over ends
00676         
00677         if(((*eiter)->GetPlexSEIdAltL()).size()!=1){
00678            MSG("UberModule",Msg::kDebug)<<"Wrong size of Plexseidaltl "
00679                                         <<((*eiter)->GetPlexSEIdAltL()).size()
00680                                         <<" snarl no "<<snarldata->snarlno
00681                                         <<" plane is "<<plane<<endl;
00682            continue;
00683         }
00684         PlexStripEndId pse = (*eiter)->GetPlexSEIdAltL().GetBestSEId();
00685         StripEnd::StripEnd_t se = pse.GetEnd();
00686         PlexSEIdAltLItem pseitem = (*eiter)->GetPlexSEIdAltL().GetBestItem();
00687 
00688 
00689         //get hit info
00690         Int_t adc = (int)(*eiter)->GetCharge();
00691         Float_t siglin = (*eiter)->GetCharge(CalDigitType::kSigLin);
00692         Float_t npe = pseitem.GetPE();
00693         Float_t sigcorr = (*eiter)->GetCharge(CalDigitType::kSigCorr);
00694         Float_t mip = Calibrator::Instance().GetMIP((*eiter)
00695                                                     ->GetCharge(CalDigitType::kSigCorr),pse);
00696         Float_t time = ((*eiter)->GetSubtractedTime(CalTimeType::kT0)
00697                         /Munits::nanosecond);
00698         PlexPixelSpotId pps=ph.GetPixelSpotId(pse);
00699         Int_t agg = pps.GetEncoded();
00700         // added by mak March 25, 2003
00701 
00702         if(iscosmic){
00703 //        snarldata->AddNextCosmicHitValues(se, adc, time, agg);
00704           continue;
00705         }
00706 
00707         snarldata->nhits++;  //increment nhits
00708         if(!gotsigcorrconv){
00709           // MAK -- 6/3/03 : Added sig fpe protection
00710           if(mip>0) snarldata->sigcorrconv = sigcorr/mip;
00711           else snarldata->sigcorrconv=-1;
00712           gotsigcorrconv=kTRUE;
00713         }
00714 
00715         snarldata->AddNextHitValues(se, adc, siglin, npe, mip, time, agg);
00716         MSG("UberModule",Msg::kDebug)<<"STRIP END "<<se<<" adc "<<adc<<" siglin "<<siglin<<" npe "<<npe<<" mip "<<mip<<endl;
00717         planemip+=mip;  //increment the mips deposited in current plane
00718         snarldata->summip[pse.GetPlane()]+=mip;  //increment mips deposited in plane array
00719         if(pse.GetPlane()%2==0){  //calculate center in even planes
00720           snarldata->mipweighcentereven+=mip*pse.GetStrip();
00721           totmipeven+=mip;
00722           if(pse.GetPlane()!=0){
00723              if(se==StripEnd::kPositive){
00724                 snarldata->totmipposeven+=mip;
00725              }
00726              else if(se==StripEnd::kNegative){
00727                 snarldata->totmipnegeven+=mip;
00728              }
00729           }
00730         }
00731         else{ //calculate center in odd planes
00732           snarldata->mipweighcenterodd+=mip*pse.GetStrip();
00733           totmipodd+=mip;
00734           if(se==StripEnd::kPositive){
00735              snarldata->totmipposodd+=mip;
00736           }
00737           else if(se==StripEnd::kNegative){
00738              snarldata->totmipnegodd+=mip;
00739           }
00740         }
00741 
00742         //look for max hit adc, npe, mip, time
00743         if(adc>snarldata->maxadc){
00744           snarldata->maxadc = adc;
00745         }
00746         if(npe>snarldata->maxnpe){
00747           snarldata->maxnpe = npe;
00748         }
00749         if(mip>snarldata->maxmip){
00750           snarldata->maxmip = mip;
00751         }
00752         if(time>snarldata->maxtime){
00753           snarldata->maxtime = time;
00754         }
00755 
00756         //add mip in all planes but plane 0
00757         if(pse.GetPlane()!=0){
00758           snarldata->totmip+=mip;
00759         }
00760 
00761         //find mips in plane 0
00762         if(pse.GetPlane()==0){
00763           snarldata->p0totmip+=mip;
00764           p0stripmip+=mip;
00765           if(mip>maxp0hitmip){
00766             maxp0hitmip = mip;
00767             snarldata->p0maxmiptstamp = time;
00768           } 
00769         }
00770 
00771         //find mips in plane 1
00772         if(pse.GetPlane()==1){
00773           snarldata->p1totmip+=mip;
00774           p1stripmip+=mip;
00775           if(mip>maxp1hitmip){
00776             maxp1hitmip = mip;
00777             snarldata->p1maxmiptstamp = time;
00778           } 
00779         }
00780       }//end loop over ends
00781 /*      
00782       if(fIsMC&&!iscosmic){
00783         PlexStripEndId mcpse = (*siter)->GetPlexSEIdAltL().GetBestSEId();
00784         mcpse.SetEnd(StripEnd::kNegative);
00785         REROOT_FLSDigit mcd = FLSDigitMap.find(mcpse)->second;
00786         MSG("UberModule",Msg::kDebug)<<"pse.GetPlane() "
00787                                   <<(*siter)
00788           ->GetPlexSEIdAltL().GetBestSEId().GetPlane()
00789                                   <<" mc plane "<<mcd.IPln()<<endl;
00790         mc->AddNextHitValues(mcd.HitBits(),mcd.SumETrue(),mcd.TPos(),mcd.RawB(),
00791                              mcd.RawA(), mcd.CorrB(), mcd.CorrA(), mcd.TDCB(),
00792                              mcd.TDCA(), mcd.AveDistTrueB(),mcd.AveDistTrueA());
00793         mc->mcenergydep += mcd.SumETrue();
00794         nmchitstrips++;
00795       }
00796 */
00797       //find max mip in an strip in plane 0
00798       if(p0stripmip>maxp0stripmip){ 
00799         maxp0stripmip = p0stripmip;
00800         snarldata->p0stripmaxmip = strip;
00801       }
00802       //find max mip in a strip in plane 1
00803       if(p1stripmip>maxp1stripmip){
00804         maxp1stripmip = p1stripmip;
00805         snarldata->p1stripmaxmip = strip;
00806       }
00807     }//end loop over strips
00808 
00809     snarldata->mipweighaveplane+=planemip*plane;  //calculate mip weigh ave plane
00810     if(planemip>maxplanemip){  //find showermax
00811       maxplanemip = planemip;
00812       snarldata->mipshowermax = maxplanemip;
00813       snarldata->showermax = plane;
00814     }
00815   }//end loop over planes
00816   // MAK -- sept 21 2003
00817   if((snarldata->totmip+snarldata->p0totmip)>0){
00818        snarldata->mipweighaveplane/=(snarldata->totmip
00819                                      +snarldata->p0totmip);
00820   }
00821   else snarldata->mipweighaveplane=-1;
00822 
00823   // MAK -- 6/3/03
00824   // added some sig fpe checking here
00825   if(totmipeven>0) snarldata->mipweighcentereven/=totmipeven;
00826   else snarldata->mipweighcentereven=-1;
00827   if(totmipodd>0) snarldata->mipweighcenterodd/=totmipodd;
00828   else snarldata->mipweighcenterodd=-1;
00829   
00830 
00831   //loop again to calculate radius of event:
00832   for(hi.Reset();hi.IsValid();hi.NextKey()){
00833     UShort_t plane = ((*hi)->GetPlexSEIdAltL()).GetPlane();
00834     if(plane>=CalDetConstants::PLANECONST){
00835       continue;
00836     }
00837     float center = 0.;
00838     if(plane%2==0){
00839       center=snarldata->mipweighcentereven;
00840     }
00841     else{
00842       center = snarldata->mipweighcenterodd;
00843     }
00844     for(CandDigitHandleItr siter(hi, kTRUE);
00845         siter.IsValid();siter.NextKey()){ //loop over strips
00846       for(CandDigitHandleItr eiter(siter,kTRUE);
00847           eiter.IsValid();eiter.Next()){ //loop over ends
00848         PlexStripEndId pse = (*eiter)->GetPlexSEIdAltL().GetBestSEId();
00849         if(pse.GetStrip()>=CalDetConstants::STRIPCONST){
00850           continue;
00851         }
00852         Float_t mip = Calibrator::Instance().GetMIP((*eiter)
00853                                                     ->GetCharge(CalDigitType::kSigCorr),pse);
00854 
00855         snarldata->mipweighrad += fabs(center-pse.GetStrip())*mip;
00856       }//end end loop
00857     }//end strip loop
00858   }//end planeloop
00859 
00860   //sum up the elements of summip
00861   for(int i=1;i<CalDetConstants::PLANECONST;i++){
00862      snarldata->summip[i]+=snarldata->summip[i-1];
00863   }
00864 
00865   // MAK -- 6/3/03 : Added sig fpe protection  
00866   if((snarldata->totmip+snarldata->p0totmip)>0.0) 
00867     snarldata->mipweighrad/=(snarldata->totmip+snarldata->p0totmip);
00868   else snarldata->mipweighrad=0;
00869 
00870   // MAK: 12 April 2005
00871   // moved call to FillNtpCalDetPID earlier in this code
00872   // rationale is that events without hits may still have pid
00873   //  this->FillNtpCalDetPID(snarldata);
00874 
00875 
00876   this -> FillNtpShower(snarldata);
00877   this -> FillNtpTrack(snarldata);
00878 //  this -> FillNtpEvent(snarldata);
00879 
00880   MSG("UberModule",Msg::kDebug)<<"giving fragment to mom"<<endl;  
00881   mom->AdoptFragment(snarldata);
00882 
00883 
00884   MSG("UberModule",Msg::kDebug)<<"done with this event"<<endl;
00885   return JobCResult(JobCResult::kPassed);
00886 }//end Reco()
00887 //________________________________________________________________________________
00888 
00889 void UberModule::FillNtpCalDetPID(UberRecord* ntprec)
00890 {
00891      MSG("UberModule",Msg::kDebug) 
00892           << "UberModuleModule::FillNtpCalDetPID" << endl;
00893 
00894      const CandCalDetPIDHandle* pidh 
00895           = dynamic_cast <const CandCalDetPIDHandle*> 
00896           (candrec -> FindCandHandle("CandCalDetPIDHandle"));
00897      if ( !pidh ) return; // all done
00898      
00899      if(pidh->NoOverlap()) ntprec->cpid.nov=1;
00900      else ntprec->cpid.nov=0;
00901 
00902      if(pidh->InCERTime()) ntprec->cpid.inct=1;
00903      else ntprec->cpid.inct=0;
00904      
00905      ntprec->cpid.pid=pidh->GetPIDType();
00906      
00907      ntprec->cpid.olchi2 = pidh->GetOLChi2();
00908      return;
00909 }
00910 
00911 void UberModule::FillNtpShower(UberRecord* ntprec)
00912 {
00913   //
00914   //  Purpose:  Private method used to fill shower portion of ntuple record.
00915   //
00916   //  Arguments: pointers to UberModuleRecord and CandRecord
00917   //  
00918   //  Return: status.
00919   // 
00920 
00921   MSG("UberModule",Msg::kDebug) << "UberModuleModule::FillNtpShower" << endl;
00922 
00923   const CandShowerListHandle *showerlisthandle 
00924    = dynamic_cast <const CandShowerListHandle*> 
00925      (candrec -> FindCandHandle("CandShowerListHandle"));
00926   if ( !showerlisthandle ) return; // all done
00927 
00928   Int_t nshower = 0;
00929   TIter showerItr(showerlisthandle->GetDaughterIterator());
00930 //  TClonesArray& showerarray = *(ntprec->shw);
00931 //2003-10-19, just store the shower with max # digits
00932   int maxdigits=0;
00933   while (CandShowerHandle* shower=dynamic_cast<CandShowerHandle*>
00934                                                       (showerItr())) {
00935     ntprec->nshw++;
00936     if(shower->GetNDigit()<=maxdigits){
00937        continue;
00938     }
00939     maxdigits = shower->GetNDigit();
00940     Int_t uid = shower->GetUidInt();
00941     fShowerUidMap.insert(std::make_pair(uid,nshower));
00942 
00943     ResetShwNStrip(shower->GetNStrip(),ntprec);
00944 
00945     // Uses new with placement to efficiently create slice ntp
00946 //    NtpSRShower* ntpshower 
00947 //         = new(showerarray[nshower++])NtpSRShower(shower->GetNStrip());
00948 
00949     nshower++;
00950     // Fill indices of associated strips in shower tree
00951     ntprec->shw.index = nshower - 1;
00952     ntprec->shw.ndigit = shower->GetNDigit();
00953     TIter showerstripItr(shower->GetDaughterIterator());
00954     Int_t nshowerstrip = 0;
00955 //    while ( CandShowerHandle *showerstrip = dynamic_cast<CandShowerHandle*>
00956     //                                                (showerstripItr()) ) {
00957     while (TObject* obj = showerstripItr()) {
00958        CandStripHandle* showerstrip = dynamic_cast<CandStripHandle*>(obj);
00959        if(showerstrip){
00960           Int_t plane = showerstrip->GetPlane();
00961           Int_t strip = showerstrip->GetStrip();
00962 //        Int_t uid = showerstrip->GetUidInt();
00963           Int_t uid = MakePlaneStripIndex(plane,strip);
00964           Int_t stripindex = fStripUidMap[uid]; // should use find
00965           ntprec->shw.AddStripAt(stripindex,nshowerstrip); // add index to strip
00966           nshowerstrip++;
00967        }
00968     }
00969 
00970     // Set range of planes included in slice
00971     ntprec->shw.plane.n = shower->GetNPlane();
00972     ntprec->shw.plane.nu = shower->GetNPlane(PlaneView::kU);
00973     ntprec->shw.plane.nv = shower->GetNPlane(PlaneView::kV);
00974     ntprec->shw.plane.beg = shower->GetBegPlane();
00975     ntprec->shw.plane.begu = shower->GetBegPlane(PlaneView::kU);
00976     ntprec->shw.plane.begv = shower->GetBegPlane(PlaneView::kV);
00977     ntprec->shw.plane.end = shower->GetEndPlane();
00978     ntprec->shw.plane.endu = shower->GetEndPlane(PlaneView::kU);
00979     ntprec->shw.plane.endv = shower->GetEndPlane(PlaneView::kV);
00980     // Set summed charge in shower
00981     ntprec->shw.ph.raw = shower->GetCharge(CalStripType::kNone);
00982     ntprec->shw.ph.siglin = shower->GetCharge(CalStripType::kSigLin);
00983     ntprec->shw.ph.sigcor = shower->GetCharge(CalStripType::kSigCorr);
00984     ntprec->shw.ph.pe = shower->GetCharge(CalStripType::kPE);
00985     ntprec->shw.ph.mip = shower->GetCharge(CalStripType::kMIP);
00986     ntprec->shw.ph.gev = shower->GetCharge(CalStripType::kGeV);
00987 
00988   }
00989 
00990   return;
00991 }//end FillNtpShower(UberRecord, candrec)
00992 
00993 //________________________________________________________________________________
00994 
00995 void UberModule::FillNtpTrack(UberRecord* ntprec)
00996 {
00997   //
00998   //  Purpose:  Private method used to fill track portion of ntuple record.
00999   //
01000   //  Arguments: pointers to UberModuleRecord and CandRecord
01001   //  
01002   //  Return: status.
01003   // 
01004 
01005 
01006   MSG("UberModule",Msg::kDebug) << "UberModule::FillNtpTrack" << endl;
01007 
01008   const Double_t cos45 = 0.70710678;  // used for converting from u,v to x,y
01009 
01010   const CandTrackListHandle *tracklisthandle 
01011      = dynamic_cast <const CandTrackListHandle*> 
01012      (candrec -> FindCandHandle("CandFitTrackListHandle"));
01013   if ( !tracklisthandle ) {
01014     tracklisthandle = dynamic_cast <const CandTrackListHandle*>
01015       (candrec -> FindCandHandle("CandTrackListHandle"));
01016   }
01017   if ( !tracklisthandle ) return;
01018 
01019   CandTrackSRListHandle* tracksrlisthandle=dynamic_cast<CandTrackSRListHandle*>
01020                       (candrec->FindCandHandle("CandTrackSRListHandle"));
01021 
01022   const VldContext& vld = *(candrec->GetVldContext());
01023 
01024   TIter trackItr(tracklisthandle->GetDaughterIterator());
01025 //  TClonesArray& trackarray = *(ntprec->trk);
01026 //2003-10-19, only store track with max # of planes
01027   int maxplanes=0;
01028 
01029   Int_t ntrack = 0;
01030   while (CandTrackHandle* track = dynamic_cast<CandTrackHandle*>(trackItr())){
01031      MSG("UberModule",Msg::kDebug) <<" max planes "<<maxplanes
01032                                    <<" NPlanes "<<track->GetNPlane()<< endl;
01033     ntprec->ntrk++;
01034     if(track->GetNPlane()<=maxplanes){
01035        continue;
01036     }
01037     maxplanes = track->GetNPlane();
01038  
01039     ResetTrkNStrip(track->GetNStrip(),ntprec);
01040 
01041     Int_t uid = track->GetUidInt();
01042     fTrackUidMap.insert(std::make_pair(uid,ntrack));
01043     // Uses new with placement to efficiently create event ntp
01044 //    NtpSRTrack* ntptrack 
01045 //             = new(trackarray[ntrack++])NtpSRTrack(track->GetNStrip());
01046 
01047     ntrack++;
01048     ntprec->trk.index = ntrack-1;
01049 
01050     CandTrackSRHandle *tracksr = dynamic_cast<CandTrackSRHandle*>(track);
01051     //CandFitTrackHandle *fittrack = dynamic_cast<CandFitTrackHandle*>(track);
01052     CandFitTrackSRHandle *fittracksr
01053                                 =dynamic_cast<CandFitTrackSRHandle*>(track);
01054 
01055     // Set range of planes included in track
01056     NtpSRTrackPlane& plane = ntprec->trk.plane;
01057     plane.n = track->GetNPlane();
01058     plane.nu = track->GetNPlane(PlaneView::kU);
01059     plane.nv = track->GetNPlane(PlaneView::kV);
01060     plane.beg = track->GetBegPlane();
01061     plane.begu = track->GetBegPlane(PlaneView::kU);
01062     plane.begv = track->GetBegPlane(PlaneView::kV);
01063     plane.end = track->GetEndPlane();
01064     plane.endu = track->GetEndPlane(PlaneView::kU);
01065     plane.endv = track->GetEndPlane(PlaneView::kV);
01066     if (tracksr) plane.ntrklike = tracksr->GetNTrackPlane();
01067     if (fittracksr) plane.ntrklike = fittracksr->GetNTrackPlane();
01068 
01069     // Set summed pulse height information for track
01070     NtpSRStripPulseHeight& ph = ntprec->trk.ph;
01071     ph.raw = track->GetCharge(CalStripType::kNone);
01072     ph.siglin = track->GetCharge(CalStripType::kSigLin);
01073     ph.sigcor = track->GetCharge(CalStripType::kSigCorr);
01074     ph.pe = track->GetCharge(CalStripType::kPE);
01075     ph.mip = track->GetCharge(CalStripType::kMIP);
01076     ph.gev = track->GetCharge(CalStripType::kGeV);
01077     
01078     ntprec->trk.ds = track->GetdS(); // distance from vertex to end
01079     ntprec->trk.range = track->GetRange(); // g/cm**2 from vertex to end
01080     // CPU to create track list
01081     if (tracksrlisthandle) ntprec->trk.cputime = tracksrlisthandle->GetCPUTime();
01082 
01083     this->FillNtpTrackVertex(&ntprec->trk,track);
01084     this->FillNtpTrackEnd(&ntprec->trk,track);
01085     this->FillNtpTrackLinearFit(&ntprec->trk,track);
01086 
01087     this->FillNtpTrackFidVtx(&ntprec->trk,track,vld);
01088     this->FillNtpTrackFidEnd(&ntprec->trk,track,vld);
01089     this->FillNtpTrackFidAll(&ntprec->trk,track,vld);
01090 
01091     this->FillNtpTrackMomentum(&ntprec->trk,track);
01092     this->FillNtpTrackFit(&ntprec->trk,track);
01093     this->FillNtpTrackTime(&ntprec->trk,track);
01094 
01095     // Loop over strips associated with track
01096     TIter trackstripItr(track->GetDaughterIterator());
01097     Int_t ntrackstrip = 0;
01098     // Fill indices of associated strips in track ntuple
01099 //    while ( CandStripHandle *trackstrip = dynamic_cast<CandStripHandle*>
01100     //                                              (trackstripItr())) {
01101     while (TObject* obj = (trackstripItr())) {
01102        CandStripHandle* trackstrip = dynamic_cast<CandStripHandle*>(obj);
01103        if(trackstrip){
01104           Int_t plane = trackstrip->GetPlane();
01105           Int_t strip = trackstrip->GetStrip();
01106           Int_t uid = MakePlaneStripIndex(plane,strip);
01107           Int_t stripindex = fStripUidMap[uid];
01108           ntprec->trk.AddStripAt(stripindex,ntrackstrip); // add index to track
01109           // Track positional information at plane associated with each strip
01110           Int_t iplane = trackstrip->GetPlane();
01111           ntprec->trk.stpu[ntrackstrip] = track->GetU(iplane);
01112           ntprec->trk.stpv[ntrackstrip] = track->GetV(iplane);
01113           ntprec->trk.stpx[ntrackstrip] = cos45*(ntprec->trk.stpu[ntrackstrip]
01114                                                -ntprec->trk.stpv[ntrackstrip]);
01115           ntprec->trk.stpy[ntrackstrip] = cos45*(ntprec->trk.stpu[ntrackstrip]
01116                                                +ntprec->trk.stpv[ntrackstrip]);
01117           ntprec->trk.stpz[ntrackstrip] = track->GetZ(iplane);
01118           // dS is travel distance from vertex
01119           ntprec->trk.stpds[ntrackstrip] = track->GetdS(iplane);
01120           
01121           // Fit track dependent quantities
01122           if ( fittracksr ) {
01123              ntprec->trk.stpfitchi2[ntrackstrip] = fittracksr->GetPlaneChi2(iplane);
01124              ntprec->trk.stpfitprechi2[ntrackstrip] 
01125                 = fittracksr->GetPlanePreChi2(iplane);
01126              ntprec->trk.stpfitqp[ntrackstrip] = fittracksr->GetPlaneQP(iplane);
01127              if ( fittracksr->GetBadFit(iplane) ) ntprec->trk.stpfit[ntrackstrip] = 0;
01128           }
01129 
01130           // Strip end dependent quantities
01131           for (UInt_t i = 0; i < 2; i++ ) {
01132              StripEnd::EStripEnd stripend = StripEnd::kNegative;
01133              if (i == 1) stripend = StripEnd::kPositive;
01134              if ( trackstrip->GetNDigit(stripend) > 0 ) {
01135                 Float_t sigmap = track->GetStripCharge(trackstrip,
01136                                                        CalStripType::kSigMapped,stripend);
01137                 Float_t mip = track->GetStripCharge(trackstrip,
01138                                                     CalStripType::kMIP,stripend);
01139                 Float_t gev = track->GetStripCharge(trackstrip,
01140                                                     CalStripType::kGeV,stripend);
01141                 ntprec->trk.SetPh(sigmap,mip,gev,ntrackstrip,i); 
01142                 ntprec->trk.SetTime(track->GetT(iplane,stripend),ntrackstrip,i);
01143              }
01144              
01145              // NJT 07/04
01146              // Use the new Calibrator to dig out the values of interest.
01147              float c0 = Calibrator::Instance().DecalStripToStrip(1.0,trackstrip->GetStripEndId(stripend));
01148              ntprec->trk.SetAttnC0( c0 ,ntrackstrip,i);
01149 
01150              float t0 = Calibrator::Instance().DecalTime(0.0,trackstrip->GetCharge(CalDigitType::kNone,stripend),trackstrip->GetStripEndId(stripend));
01151              ntprec->trk.SetCalT0(t0,ntrackstrip,i);
01152           }
01153           ntrackstrip++;
01154        }
01155     }
01156         
01157   }
01158 
01159   return;
01160 
01161 }//end FillNtpTrack()
01162 
01163 //________________________________________________________________________________
01164 /*
01165 void UberModule::FillNtpEvent(UberRecord* ntprec)
01166 {
01167   //
01168   //  Purpose:  Private method used to fill event portion of ntuple record.
01169   //
01170   //  Arguments: pointers to UberModuleRecord and CandRecord
01171   //  
01172   //  Return: status.
01173   // 
01174 
01175   MSG("UberModule",Msg::kDebug) << "UberModule::FillNtpEvent" << endl;
01176 
01177   const CandEventListHandle *eventlisthandle 
01178    = dynamic_cast <const CandEventListHandle*> 
01179      (candrec -> FindCandHandle("CandEventListHandle"));
01180   if ( !eventlisthandle ) return; // all done
01181 
01182   Int_t nevent = 0;
01183   TIter eventItr(eventlisthandle->GetDaughterIterator());
01184 //  TClonesArray& eventarray = *(ntprec->evt);
01185 //2003-10-19, only store track with max # of digits
01186   int maxdigits=0;
01187 
01188   while (CandEventHandle* event=dynamic_cast<CandEventHandle*>
01189                                                       (eventItr())) {
01190      ntprec->nevt++;
01191     if(event->GetNDigit()<=maxdigits){
01192        continue;
01193     }
01194     maxdigits = event->GetNDigit();
01195     if(ntprec->evt){
01196        ntprec->evt->Clear();
01197        delete ntprec->evt;
01198        ntprec->evt = 0;
01199     }
01200     ntprec->evt = new NtpSREvent(event->GetNStrip(),
01201                                  event->GetLastShower()+1,
01202                                  event->GetLastTrack()+1);
01203 
01204     Int_t uid = event->GetUidInt();
01205     fEventUidMap.insert(std::make_pair(uid,nevent));
01206     // Uses new with placement to efficiently create event ntp
01207 //    NtpSREvent* ntpevent 
01208 //    = new(eventarray[nevent++])NtpSREvent(event->GetNStrip(),
01209 //                      event->GetLastShower()+1,event->GetLastTrack()+1);
01210 
01211     ntprec->evt->index = nevent-1;
01212     ntprec->evt->ndigit = event->GetNDigit();
01213 
01214     TIter eventstripItr(event->GetDaughterIterator());
01215     Int_t neventstrip = 0;
01216     // Fill indices of associated strips,showers,tracks in event ntuple
01217     while ( CandStripHandle *eventstrip = dynamic_cast<CandStripHandle*>
01218                                                       (eventstripItr())) {
01219       Int_t stripindex = fStripUidMap[eventstrip->GetUidInt()];
01220       ntprec->evt->AddStripAt(stripindex,neventstrip); // add index to event
01221       neventstrip++;
01222     }
01223     for (Int_t i = 0; i <= event->GetLastTrack(); i++ ) {
01224       const CandTrackHandle* track = event->GetTrack(i);
01225       Int_t trackindex = fTrackUidMap[track->GetUidInt()];
01226       ntprec->evt -> AddTrackAt(trackindex,i);
01227     }
01228     for (Int_t i = 0; i <= event->GetLastShower(); i++ ) {
01229       const CandShowerHandle* shower = event->GetShower(i);
01230       Int_t showerindex = fShowerUidMap[shower->GetUidInt()];
01231       ntprec->evt -> AddShowerAt(showerindex,i);
01232     }
01233 
01234     // Set range of planes included in event
01235     ntprec->evt->plane.n = event->GetNPlane();
01236     ntprec->evt->plane.nu = event->GetNPlane(PlaneView::kU);
01237     ntprec->evt->plane.nv = event->GetNPlane(PlaneView::kV);
01238     ntprec->evt->plane.beg = event->GetBegPlane();
01239     ntprec->evt->plane.begu = event->GetBegPlane(PlaneView::kU);
01240     ntprec->evt->plane.begv = event->GetBegPlane(PlaneView::kV);
01241     ntprec->evt->plane.end = event->GetEndPlane();
01242     ntprec->evt->plane.endu = event->GetEndPlane(PlaneView::kU);
01243     ntprec->evt->plane.endv = event->GetEndPlane(PlaneView::kV);
01244     // Set summed charge in event
01245     ntprec->evt->ph.raw = event->GetCharge(CalStripType::kNone);
01246     ntprec->evt->ph.siglin = event->GetCharge(CalStripType::kSigLin);
01247     ntprec->evt->ph.sigcor = event->GetCharge(CalStripType::kSigCorr);
01248     ntprec->evt->ph.pe = event->GetCharge(CalStripType::kPE);
01249     ntprec->evt->ph.mip = event->GetCharge(CalStripType::kMIP);
01250     ntprec->evt->ph.gev = event->GetCharge(CalStripType::kGeV);
01251   }
01252 
01253   return;
01254 }//end FillNtpEvent()
01255 */
01256 //________________________________________________________________________________
01257 
01258 
01259 void UberModule::FillNtpTrackMomentum(NtpSRTrack* ntptrack,
01260                                        const CandTrackHandle* track) {
01261   //
01262   //  Purpose:  Private method used to fill NtpSRTrack momentum data member
01263   //            given a pointer to the ntuple track and a pointer to the
01264   //            associated CandTrackHandle.
01265   //
01266   //  Return: none.
01267   // 
01268 
01269 
01270   MSG("UberModule",Msg::kDebug) << "UberModule::FillNtpTrackMomentum" << endl;
01271 
01272   NtpSRMomentum& momentum = ntptrack->momentum;
01273   momentum.range = track->GetMomentum();
01274 
01275   const CandFitTrackHandle* fittrack 
01276     = dynamic_cast<const CandFitTrackHandle*>(track);
01277   if ( fittrack ) {
01278       if(fittrack->GetMomentumCurve())
01279           momentum.qp = fittrack->GetEMCharge()/fittrack->GetMomentumCurve();
01280       else momentum.qp=0;
01281   }
01282 
01283   const CandFitTrackSRHandle* fittracksr 
01284     = dynamic_cast<const CandFitTrackSRHandle*>(track);
01285   if ( fittracksr ) {
01286     momentum.eqp = fittracksr->GetVtxQPError();
01287   }
01288 
01289   return;
01290 
01291 }
01292 
01293 void UberModule::FillNtpTrackTime(NtpSRTrack* ntptrack,
01294                                    const CandTrackHandle* track) {
01295   //
01296   //  Purpose:  Private method used to fill NtpSRTrack time data member
01297   //            given a pointer to the ntuple track and a pointer to the
01298   //            associated CandTrackHandle.
01299   //
01300   //  Return: none.
01301   // 
01302 
01303 
01304   MSG("UberModule",Msg::kDebug) << "UberModule::FillNtpTrackTime" << endl;
01305 
01306   NtpSRTrackTime& time = ntptrack->time;
01307 
01308   const CandTrackSRHandle* tracksr 
01309     = dynamic_cast<const CandTrackSRHandle*>(track);
01310   if ( tracksr ) {
01311     time.ndigit = tracksr->GetNTimeFitDigit();
01312     time.chi2   = tracksr->GetTimeFitChi2();
01313   }
01314 
01315   const CandFitTrackSRHandle* fittracksr 
01316     = dynamic_cast<const CandFitTrackSRHandle*>(track);
01317   if ( fittracksr ) {
01318     time.ndigit = fittracksr->GetNTimeFitDigit();
01319     time.chi2   = fittracksr->GetTimeFitChi2();
01320   }
01321 
01322   time.cdtds = fabs(track->GetTimeSlope())*3.e8;
01323   time.dtds = track->GetTimeSlope();
01324   time.t0   = track->GetTimeOffset();
01325   
01326   Double_t totuph[2] = {0}, totvph[2] = {0};
01327   Int_t    ndutimespace[1000] = {0}, ndvtimespace[1000] = {0};
01328   Double_t dutimespace[1000]  = {0}, dvtimespace[1000] = {0};
01329   
01330   TIter trackstripItr(track->GetDaughterIterator());
01331   while (CandStripHandle* trackstrip = dynamic_cast<CandStripHandle*>
01332                                                    (trackstripItr())) {
01333     Int_t iplane = trackstrip -> GetPlane();
01334     assert(iplane >= 0 && iplane < 1000);
01335     Float_t ph0 = trackstrip -> GetCharge(StripEnd::kNegative);
01336     Float_t ph1 = trackstrip -> GetCharge(StripEnd::kPositive);
01337     Float_t t0  = track->GetT(iplane,StripEnd::kNegative);
01338     Float_t t1  = track->GetT(iplane,StripEnd::kPositive);
01339 
01340     if ( trackstrip -> GetPlaneView() == PlaneView::kU ) {
01341       if ( trackstrip -> GetNDigit(StripEnd::kNegative) > 0 ) {
01342         totuph[0] += ph0;
01343         time.u0   += ph0*t0;
01344       }
01345       if ( trackstrip -> GetNDigit(StripEnd::kPositive) > 0 ) {
01346         totuph[1] += ph1;
01347         time.u1   += ph1*t1;
01348       }
01349       if ( t0 > 0 && t1 > 0 ) {
01350         Float_t dapos = .5*(t0-t1)*PropagationVelocity::VelocityPosErr();
01351         if (!ndutimespace[iplane] || fabs(dapos)<fabs(dutimespace[iplane])) {
01352           dutimespace[iplane] = dapos;
01353           ndutimespace[iplane] = 1;
01354         }
01355       }
01356     }
01357 
01358     if ( trackstrip -> GetPlaneView() == PlaneView::kV ) {
01359       if ( trackstrip -> GetNDigit(StripEnd::kNegative) > 0 ) {
01360         totvph[0] += ph0;
01361         time.v0   += ph0*t0;
01362       }
01363       if ( trackstrip -> GetNDigit(StripEnd::kPositive) > 0 ) {
01364         totvph[1] += ph1;
01365         time.v1   += ph1*t1;
01366       }
01367       if ( trackstrip -> GetNDigit(StripEnd::kNegative) > 0 &&
01368            trackstrip -> GetNDigit(StripEnd::kPositive) > 0 ) {
01369         Float_t dapos = .5*(t0-t1)*PropagationVelocity::VelocityPosErr();
01370         if (!ndvtimespace[iplane]||fabs(dapos) < fabs(dvtimespace[iplane])) {
01371           dvtimespace[iplane] = dapos;
01372           ndvtimespace[iplane] = 1;
01373         }
01374       }
01375     }
01376   }
01377 
01378   Int_t ndu = 0;
01379   Int_t ndv = 0;
01380   for (int ip = 0; ip < 1000; ip++ ) {
01381     if ( ndutimespace[ip] > 0 ) {
01382       ndu++;
01383       time.du += dutimespace[ip];
01384     }
01385     if ( ndvtimespace[ip] > 0 ) {
01386       ndv++;
01387       time.dv += dvtimespace[ip];
01388     }
01389   }
01390   if ( totuph[0] > 0. ) time.u0 /= totuph[0];
01391   if ( totvph[0] > 0. ) time.v0 /= totvph[0];
01392   if ( totuph[1] > 0. ) time.u1 /= totuph[1];
01393   if ( totvph[1] > 0. ) time.v1 /= totvph[1];
01394   if ( ndu ) time.du /= (Float_t)ndu;
01395   if ( ndv ) time.dv /= (Float_t)ndv;
01396 
01397   return;
01398 
01399 }
01400 
01401 void UberModule::FillNtpTrackFit(NtpSRTrack* ntptrack,
01402                                   const CandTrackHandle* track) {
01403   //
01404   //  Purpose:  Private method used to fill NtpSRTrack fit data member
01405   //            given a pointer to the ntuple track and a pointer to the
01406   //            associated CandTrackHandle.
01407   //
01408   //  Return: none.
01409   // 
01410 
01411 
01412   MSG("UberModule",Msg::kDebug) << "UberModule::FillNtpTrackFit" << endl;
01413 
01414   NtpSRFitTrack& fit = ntptrack->fit;
01415 
01416   const CandFitTrackHandle* fittrack 
01417     = dynamic_cast<const CandFitTrackHandle*>(track);
01418   if ( fittrack ) {
01419     fit.chi2 = fittrack->GetChi2();
01420     fit.pass = fittrack->GetPass();
01421   }
01422 
01423   const CandFitTrackSRHandle* fittracksr 
01424     = dynamic_cast<const CandFitTrackSRHandle*>(track);
01425   if ( fittracksr ) {
01426     fit.ndof      = fittracksr->GetNDOF();
01427     fit.niterate  = fittracksr->GetNIterate();
01428     fit.nswimfail = fittracksr->GetNSwimFail();
01429     fit.cputime   = fittracksr->GetCPUTime();
01430   }
01431 
01432   return;
01433 
01434 }
01435 
01436 void UberModule::FillNtpTrackFidVtx(NtpSRTrack* ntptrack,
01437                                      const CandTrackHandle* track,
01438                                      const VldContext& vld) {
01439   //
01440   //  Purpose:  Private method used to fill NtpSRTrack fidvtx data member
01441   //            given a pointer to the ntuple track, a pointer to the
01442   //            associated CandTrackHandle, and the vldcontext of the
01443   //            event.
01444   //
01445   //  Return: none.
01446   // 
01447 
01448 
01449   MSG("UberModule",Msg::kDebug) << "UberModule::FillNtpTrackFidVtx" << endl;
01450 
01451   NtpSRFiducial& fid = ntptrack->fidvtx;
01452   NtpSRVertex& vtx = ntptrack->vtx;
01453   //NtpSRFitVertex& vtx = ntptrack->vtx;
01454   this->FillNtpFiducial(fid,vtx,vld); // fills dr,dz
01455 
01456   fid.trace   = track->GetVtxTrace();
01457   fid.tracez  = track->GetVtxTraceZ();
01458 
01459   const CandFitTrackSRHandle* fittracksr 
01460     = dynamic_cast<const CandFitTrackSRHandle*>(track);
01461   if ( fittracksr ) {                                     
01462     fid.nplane  = fittracksr->GetVtxExtrapolate();
01463     //   fid.nplaneu = fittracksr->GetVtxExtrapolate(PlaneView::kU);
01464     //  fid.nplanev = fittracksr->GetVtxExtrapolate(PlaneView::kV);
01465   }  
01466 
01467   return;
01468 
01469 }
01470 
01471 void UberModule::FillNtpTrackFidEnd(NtpSRTrack* ntptrack,
01472                                      const CandTrackHandle* track,
01473                                      const VldContext& vld) {
01474   //
01475   //  Purpose:  Private method used to fill NtpSRTrack fidend data member
01476   //            given a pointer to the ntuple track, a pointer to the
01477   //            associated CandTrackHandle, and the vldcontext of the
01478   //            event.
01479   //
01480   //  Return: none.
01481   // 
01482 
01483 
01484   MSG("UberModule",Msg::kDebug) << "UberModule::FillNtpTrackFidEnd" << endl;
01485 
01486   NtpSRFiducial& fid = ntptrack->fidend;
01487   NtpSRVertex& end = ntptrack->end;
01488   //NtpSRFitVertex& end = ntptrack->end;
01489   this->FillNtpFiducial(fid,end,vld); // fills dr,dz
01490 
01491   fid.trace   = track->GetEndTrace();
01492   fid.tracez  = track->GetEndTraceZ();
01493 
01494   const CandFitTrackSRHandle* fittracksr 
01495     = dynamic_cast<const CandFitTrackSRHandle*>(track);
01496   if ( fittracksr ) {                                     
01497     fid.nplane  = fittracksr->GetEndExtrapolate();
01498     //  fid.nplaneu = fittracksr->GetEndExtrapolate(PlaneView::kU);
01499     //  fid.nplanev = fittracksr->GetEndExtrapolate(PlaneView::kV);
01500   }  
01501 
01502   return;
01503 
01504 }
01505 
01506 void UberModule::FillNtpTrackFidAll(NtpSRTrack* ntptrack,
01507                                      const CandTrackHandle* track,
01508                                      const VldContext& vld) {
01509   //
01510   //  Purpose:  Private method used to fill NtpSRTrack fidall data member
01511   //            given a pointer to the ntuple track, a pointer to the
01512   //            associated CandTrackHandle, and the vldcontext of the
01513   //            event.
01514   //
01515   //  Return: none.
01516   //
01517   //  Notes: This method should be called after FillNtpTrackFidAll
01518   //         and FillNtpFiducialEnd 
01519 
01520 
01521   MSG("UberModule",Msg::kDebug) << "UberModule::FillNtpFiducialAll" << endl;
01522   const Double_t cos45 = 0.70710678;  // used for converting from u,v to x,y
01523 
01524   NtpSRFiducial& fidall = ntptrack->fidall;
01525 
01526   NtpSRFiducial fidtmp;
01527   NtpSRVertex vtxtmp;
01528   fidall.dr = min(ntptrack->fidvtx.dr,ntptrack->fidend.dr);
01529   fidall.dz = min(ntptrack->fidvtx.dz,ntptrack->fidend.dz);
01530   fidall.trace   = min(ntptrack->fidvtx.trace,ntptrack->fidend.trace);
01531   fidall.tracez  = min(ntptrack->fidvtx.tracez,ntptrack->fidend.tracez);
01532   fidall.nplane  = min(ntptrack->fidvtx.nplane,ntptrack->fidend.nplane);
01533   // fidall.nplaneu  = min(ntptrack->fidvtx.nplaneu,ntptrack->fidend.nplaneu);
01534   // fidall.nplanev  = min(ntptrack->fidvtx.nplanev,ntptrack->fidend.nplanev);
01535 
01536   TIter trackstripItr(track->GetDaughterIterator());
01537   while(CandStripHandle* trackstrip = dynamic_cast<CandStripHandle*>
01538         (trackstripItr())) {
01539     // For all strips calculate closest approach to boundaries
01540     Int_t iplane = trackstrip->GetPlane();
01541     vtxtmp.plane = iplane;
01542     vtxtmp.u = track->GetU(iplane);
01543     vtxtmp.v = track->GetV(iplane);
01544     vtxtmp.x = cos45*(vtxtmp.u-vtxtmp.v);
01545     vtxtmp.y = cos45*(vtxtmp.u+vtxtmp.v);
01546     this->FillNtpFiducial(fidtmp,vtxtmp,vld); // fills dr,dz
01547     fidall.dr = min(fidtmp.dr,fidall.dr);
01548     fidall.dz = min(fidtmp.dz,fidall.dz);
01549   }
01550 
01551   return;
01552 
01553 }
01554 
01555 void UberModule::FillNtpTrackVertex(NtpSRTrack* ntptrack,
01556                                      const CandTrackHandle* track) {
01557   //
01558   //  Purpose:  Private method used to fill vertex portion of ntuple track.
01559   //
01560   //  Arguments: NtpSRTrack to fill and const CandTrackHandle pointer.
01561   //  
01562   //  Return: status.
01563   // 
01564 
01565 
01566   MSG("UberModule",Msg::kDebug) << "UberModule::FillNtpTrackVertex" << endl;
01567 
01568   const Double_t cos45 = 0.70710678;  // used for converting from u,v to x,y
01569 
01570   NtpSRVertex& vtx = ntptrack->vtx;
01571   //  NtpSRFitVertex& vtx = ntptrack->vtx;
01572 
01573   vtx.u = track->GetVtxU();
01574   vtx.v = track->GetVtxV();
01575   vtx.x = cos45*(vtx.u-vtx.v);
01576   vtx.y = cos45*(vtx.u+vtx.v);
01577   vtx.z = track->GetVtxZ();
01578   vtx.t = track->GetVtxT();
01579   vtx.plane = track->GetVtxPlane();
01580   vtx.dcosu = track->GetDirCosU();
01581   vtx.dcosv = track->GetDirCosV();
01582   vtx.dcosx = cos45*(vtx.dcosu-vtx.dcosv);
01583   vtx.dcosy = cos45*(vtx.dcosu+vtx.dcosv);
01584   vtx.dcosz = track->GetDirCosZ();
01585 
01586   const CandFitTrackSRHandle* fittracksr
01587                         = dynamic_cast<const CandFitTrackSRHandle*>(track);
01588   if ( fittracksr ) {
01589     vtx.eu = fittracksr->GetVtxUError();
01590     vtx.ev = fittracksr->GetVtxVError();
01591     vtx.ex = cos45*sqrt(vtx.eu*vtx.eu+vtx.ev*vtx.ev);
01592     vtx.ey = cos45*sqrt(vtx.eu*vtx.eu+vtx.ev*vtx.ev);
01593     Double_t edudz = fittracksr->GetVtxdUError();
01594     Double_t edvdz = fittracksr->GetVtxdVError();
01595     // These calculations should include the dudz and dvdz covariance terms
01596     // but currently the covariance terms are not accessible
01597     vtx.edcosz = fabs(vtx.dcosz)*sqrt(pow(vtx.dcosu*edudz,2)+
01598                                       pow(vtx.dcosv*edvdz,2));
01599     vtx.edcosu = sqrt(fabs(vtx.dcosz))*sqrt(pow(vtx.dcosu*vtx.dcosv*edvdz,2)
01600                + pow((pow(vtx.dcosz,2)+pow(vtx.dcosv,2))*edudz,2));
01601     vtx.edcosv = sqrt(fabs(vtx.dcosz))*sqrt(pow(vtx.dcosu*vtx.dcosv*edudz,2)
01602                + pow((pow(vtx.dcosz,2)+pow(vtx.dcosu,2))*edvdz,2));
01603     vtx.edcosx = cos45*sqrt(vtx.edcosu*vtx.edcosu+vtx.edcosv*vtx.edcosv);
01604     vtx.edcosy = cos45*sqrt(vtx.edcosu*vtx.edcosu+vtx.edcosv*vtx.edcosv);
01605   }
01606 
01607   return;
01608 
01609 }
01610 
01611 void UberModule::FillNtpTrackEnd(NtpSRTrack* ntptrack,
01612                                   const CandTrackHandle* track) {
01613   //
01614   //  Purpose:  Private method used to fill end portion of ntuple track.
01615   //
01616   //  Arguments: NtpSRTrack to fill and const CandTrackHandle pointer.
01617   //  
01618   //  Return: status.
01619   // 
01620 
01621 
01622   MSG("UberModule",Msg::kDebug) << "UberModule::FillNtpTrackEnd" << endl;
01623 
01624   const Double_t cos45 = 0.70710678;  // used for converting from u,v to x,y
01625 
01626   NtpSRVertex& end = ntptrack->end;
01627   //  NtpSRFitVertex& end = ntptrack->end;
01628 
01629   end.u = track->GetEndU();
01630   end.v = track->GetEndV();
01631   end.x = cos45*(end.u-end.v);
01632   end.y = cos45*(end.u+end.v);
01633   end.z = track->GetEndZ();
01634   end.t = track->GetEndT();
01635   end.plane = track->GetEndPlane();
01636   end.dcosu = track->GetEndDirCosU();
01637   end.dcosv = track->GetEndDirCosV();
01638   end.dcosx = cos45*(end.dcosu-end.dcosv);
01639   end.dcosy = cos45*(end.dcosu+end.dcosv);
01640   end.dcosz = track->GetEndDirCosZ();
01641 
01642   const CandFitTrackSRHandle* fittracksr
01643                     =dynamic_cast<const CandFitTrackSRHandle*>(track);
01644   if ( fittracksr ) {
01645     end.eu = fittracksr->GetEndUError();
01646     end.ev = fittracksr->GetEndVError();
01647     end.ex = cos45*sqrt(end.eu*end.eu+end.ev*end.ev);
01648     end.ey = cos45*sqrt(end.eu*end.eu+end.ev*end.ev);
01649     Double_t edudz = fittracksr->GetEnddUError();
01650     Double_t edvdz = fittracksr->GetEnddVError();
01651     // These calculations should include the dudz and dvdz covariance terms
01652     // but currently the covariance terms are not accessible
01653     end.edcosz = fabs(end.dcosz)*sqrt(pow(end.dcosu*edudz,2)+
01654                                       pow(end.dcosv*edvdz,2));
01655     end.edcosu = sqrt(fabs(end.dcosz))*sqrt(pow(end.dcosu*end.dcosv*edvdz,2)
01656                + pow((pow(end.dcosz,2)+pow(end.dcosv,2))*edudz,2));
01657     end.edcosv = sqrt(fabs(end.dcosz))*sqrt(pow(end.dcosu*end.dcosv*edudz,2)
01658                + pow((pow(end.dcosz,2)+pow(end.dcosu,2))*edvdz,2));
01659     end.edcosx = cos45*sqrt(end.edcosu*end.edcosu+end.edcosv*end.edcosv);
01660     end.edcosy = cos45*sqrt(end.edcosu*end.edcosu+end.edcosv*end.edcosv);
01661   }
01662 
01663   return;
01664 
01665 }
01666 
01667 void UberModule::FillNtpTrackLinearFit(NtpSRTrack* ntptrack,
01668                                         const CandTrackHandle* track) {
01669   //
01670   //  Purpose:  Private method used to fill linar fit track portion of 
01671   //            ntuple track.
01672   //
01673   //  Arguments: NtpSRTrack to fill and const CandTrackHandle pointer.
01674   //  
01675   //  Return: status.
01676   // 
01677 
01678 
01679   MSG("UberModule",Msg::kDebug) << "UberModule::FillNtpFitTrack" << endl;
01680 
01681   const Double_t cos45 = 0.70710678;  // used for converting from u,v to x,y
01682 
01683   // First array dimension is view (u or v)
01684   // Second array dimension is plane number
01685   const Int_t nplane = 1000;
01686   const Int_t nview = 2; // u,v
01687   Double_t zfit[nview][nplane]={{0},{0}},fit[nview][nplane]={{0},{0}};
01688   Double_t wfit[nview][nplane]={{0},{0}},ph[nview][nplane]={{0},{0}};
01689 
01690   // Loop over all strips on track
01691   TIter trackstripItr(track->GetDaughterIterator());
01692   while(CandStripHandle* trackstrip = dynamic_cast<CandStripHandle*>
01693         (trackstripItr())) {
01694 
01695     Float_t tpos = trackstrip->GetTPos();
01696     Float_t phend[2]  = {trackstrip->GetCharge(StripEnd::kNegative),
01697                          trackstrip->GetCharge(StripEnd::kPositive)};
01698     Float_t phsum   = phend[0] + phend[1];
01699 
01700     Int_t iplane = trackstrip->GetPlane();
01701     assert(iplane >= 0 && iplane < 1000);
01702 
01703     Int_t iview = -1;
01704     if ( trackstrip->GetPlaneView()==PlaneView::kU ) iview = 0;
01705     else if ( trackstrip->GetPlaneView()==PlaneView::kV ) iview = 1;
01706     else continue;
01707 
01708     zfit[iview][iplane] = trackstrip->GetZPos();
01709     fit[iview][iplane] += tpos*phsum; // pulse height weighted average of tpos
01710     ph[iview][iplane]  += phsum;  // summed pulse height both ends 
01711     wfit[iview][iplane] = 1;   // weight assigned for use in fit
01712   }
01713 
01714   // Finished with loop over all strips, calculate fit results
01715   for ( Int_t ip = 0; ip < nplane; ip++ ) {
01716     for ( Int_t iv = 0; iv < nview; iv++ ) {
01717       if ( ph[iv][ip] > 0 ) fit[iv][ip] /= ph[iv][ip];
01718     }
01719   }
01720 
01721   Double_t uparm[2],ueparm[2],vparm[2],veparm[2];
01722   LinearFit::Weighted(1000,zfit[0],fit[0],wfit[0],uparm,ueparm);
01723   LinearFit::Weighted(1000,zfit[1],fit[1],wfit[1],vparm,veparm);
01724   NtpSRVertex& lin = ntptrack->lin;
01725 
01726   Double_t dudz = uparm[1];
01727   Double_t dvdz = vparm[1];
01728   Double_t dxdz = cos45*(uparm[1]-vparm[1]);
01729   Double_t dydz = cos45*(uparm[1]+vparm[1]);
01730   lin.z = track->GetVtxZ();
01731   lin.t = track->GetVtxT();
01732   lin.u = uparm[0]+lin.z*dudz;
01733   lin.v = vparm[0]+lin.z*dvdz;
01734   lin.x = cos45*(uparm[0]-vparm[0])+lin.z*dxdz;
01735   lin.y = cos45*(uparm[0]+vparm[0])+lin.z*dydz;
01736   Double_t anorm = sqrt(1.+dxdz*dxdz+dydz*dydz);
01737   Int_t sign = 1;
01738   if ( track->GetTimeSlope() < 0. ) sign = -1;
01739   lin.dcosu = sign*dudz/anorm;
01740   lin.dcosv = sign*dvdz/anorm;
01741   lin.dcosx = sign*dxdz/anorm;
01742   lin.dcosy = sign*dydz/anorm;
01743   lin.dcosz = sign*1./anorm;
01744 
01745   return;
01746 
01747 }
01748 
01749 
01750 
01751 
01752 
01753 void UberModule::FillNtpFiducial(NtpSRFiducial& fid, const NtpSRVertex& vtx, 
01754                                   const VldContext& vld) {
01755   //
01756   //  Purpose:  Private method used to fill NtpSRFiducial given a point
01757   //            defined by NtpSRVertex and the detector type.
01758   //
01759   //  Return: none.
01760   // 
01761 
01762 
01763   MSG("UberModule",Msg::kDebug) << "UberModule::FillNtpFiducial" << endl;
01764 
01765   UgliGeomHandle ugh(vld);
01766   Float_t zextent[2];
01767   ugh.GetZExtent(zextent[0],zextent[1]); // apparently this doesn't work yet?
01768   Detector::Detector_t dtype = vld.GetDetector();
01769 
01770   Float_t du,dv,dx,dy;
01771 
01772   switch (dtype) {
01773 
01774   case Detector::kFar:
01775     du = min(4.-vtx.u,4.+vtx.u);
01776     dv = min(4.-vtx.v,4.+vtx.v);
01777     dx = min(4.-vtx.x,4.+vtx.x);
01778     dy = min(4.-vtx.y,4.+vtx.y);
01779     fid.dr = min(min(du,dv),min(dx,dy));
01780     if ( fid.dr < 0. ) fid.dr = 0;
01781     fid.dz = min(vtx.z-zextent[0],zextent[1]-vtx.z);
01782     
01783     break;
01784 
01785   case Detector::kCalDet:
01786     fid.dr = min(min(0.5+vtx.x,0.5-vtx.x),min(0.5+vtx.y,0.5-vtx.y));
01787     if ( fid.dr < 0. ) fid.dr = 0;
01788     fid.dz = min(vtx.z-zextent[0],zextent[1]-vtx.z);
01789 
01790     break;
01791 
01792   default:
01793     MSG("EventSR",Msg::kWarning) << "Detector type " << dtype
01794                                  << " not supported." << endl;
01795     break;
01796   }
01797 
01798   return;
01799 
01800 }
01801 //_____________________________________________________________________________
01802 void UberModule::ResetShwNStrip(Int_t nstrip, UberRecord *ntprec)
01803 {
01804   // MAK: Feb 15, 2005
01805   // rewrote this routine, so it doesn't delete and call new every event
01806 
01807   Int_t old_nstrip = ntprec->shw.nstrip;
01808   if(old_nstrip<nstrip){ 
01809     ntprec->shw.Clear();// deletes shw.stp
01810     ntprec->shw.stp = new Int_t[nstrip];
01811     for ( Int_t i = 0; i < nstrip; i++ ) ntprec->shw.stp[i] = -1;
01812   }
01813   else{
01814     for ( Int_t i = 0; i < old_nstrip; i++ ) ntprec->shw.stp[i] = -1;
01815   }
01816   ntprec->shw.nstrip = nstrip;
01817 
01818 }
01819 //_____________________________________________________________________________
01820 void UberModule::ResetTrkNStrip(Int_t nstrip, UberRecord *ntprec)
01821 {
01822   // MAK: Feb 15, 2005
01823   // rewrote this routine, so it doesn't delete and call new every event
01824   // and so it DOES call delete before calling new
01825 
01826   Int_t old_nstrip = ntprec->trk.nstrip;
01827   if ( old_nstrip < nstrip ) {
01828     ntprec->trk.Clear();
01829     ntprec->trk.stp           = new Int_t[nstrip]; 
01830     ntprec->trk.stpfit        = new Byte_t[nstrip]; 
01831     ntprec->trk.stpu          = new Float_t[nstrip]; 
01832     ntprec->trk.stpv          = new Float_t[nstrip]; 
01833     ntprec->trk.stpx          = new Float_t[nstrip]; 
01834     ntprec->trk.stpy          = new Float_t[nstrip]; 
01835     ntprec->trk.stpz          = new Float_t[nstrip]; 
01836     ntprec->trk.stpds         = new Float_t[nstrip]; 
01837     ntprec->trk.stpfitchi2    = new Float_t[nstrip]; 
01838     ntprec->trk.stpfitprechi2 = new Float_t[nstrip]; 
01839     ntprec->trk.stpfitqp      = new Float_t[nstrip]; 
01840     ntprec->trk.stpph0sigmap  = new Float_t[nstrip]; 
01841     ntprec->trk.stpph0mip     = new Float_t[nstrip]; 
01842     ntprec->trk.stpph0gev     = new Float_t[nstrip]; 
01843     ntprec->trk.stpph1sigmap  = new Float_t[nstrip]; 
01844     ntprec->trk.stpph1mip     = new Float_t[nstrip]; 
01845     ntprec->trk.stpph1gev     = new Float_t[nstrip]; 
01846     ntprec->trk.stpattn0c0    = new Float_t[nstrip]; 
01847     ntprec->trk.stpattn1c0    = new Float_t[nstrip]; 
01848     ntprec->trk.stpt0         = new Double_t[nstrip]; 
01849     ntprec->trk.stpt1         = new Double_t[nstrip]; 
01850     ntprec->trk.stptcal0t0    = new Double_t[nstrip]; 
01851     ntprec->trk.stptcal1t0    = new Double_t[nstrip]; 
01852 
01853     for (Int_t i = 0; i < nstrip; i++ ) {
01854       ntprec->trk.stp[i] = -1;
01855       ntprec->trk.stpfit[i] = 1;
01856       ntprec->trk.stpu[i] = -999999;
01857       ntprec->trk.stpv[i] = -999999;
01858       ntprec->trk.stpx[i] = -999999;
01859       ntprec->trk.stpy[i] = -999999;
01860       ntprec->trk.stpz[i] = -999999;
01861       ntprec->trk.stpds[i] = -999999;
01862       ntprec->trk.stpfitchi2[i] = 0;
01863       ntprec->trk.stpfitprechi2[i] = 0;
01864       ntprec->trk.stpfitqp[i] = 0;
01865       ntprec->trk.stpph0sigmap[i] = -999999;
01866       ntprec->trk.stpph0mip[i] = -999999;
01867       ntprec->trk.stpph0gev[i] = -999999;
01868       ntprec->trk.stpph1sigmap[i] = -999999;
01869       ntprec->trk.stpph1mip[i] = -999999;
01870       ntprec->trk.stpph1gev[i] = -999999;
01871       ntprec->trk.stpattn0c0[i] = -999999;
01872       ntprec->trk.stpattn1c0[i] = -999999;
01873       ntprec->trk.stpt0[i] = -999999;
01874       ntprec->trk.stpt1[i] = -999999;
01875       ntprec->trk.stptcal0t0[i] = -999999;
01876       ntprec->trk.stptcal1t0[i] = -999999;
01877     }
01878   }
01879   else{
01880     for (Int_t i = 0; i < old_nstrip; i++ ) {
01881       ntprec->trk.stp[i] = -1;
01882       ntprec->trk.stpfit[i] = 1;
01883       ntprec->trk.stpu[i] = -999999;
01884       ntprec->trk.stpv[i] = -999999;
01885       ntprec->trk.stpx[i] = -999999;
01886       ntprec->trk.stpy[i] = -999999;
01887       ntprec->trk.stpz[i] = -999999;
01888       ntprec->trk.stpds[i] = -999999;
01889       ntprec->trk.stpfitchi2[i] = 0;
01890       ntprec->trk.stpfitprechi2[i] = 0;
01891       ntprec->trk.stpfitqp[i] = 0;
01892       ntprec->trk.stpph0sigmap[i] = -999999;
01893       ntprec->trk.stpph0mip[i] = -999999;
01894       ntprec->trk.stpph0gev[i] = -999999;
01895       ntprec->trk.stpph1sigmap[i] = -999999;
01896       ntprec->trk.stpph1mip[i] = -999999;
01897       ntprec->trk.stpph1gev[i] = -999999;
01898       ntprec->trk.stpattn0c0[i] = -999999;
01899       ntprec->trk.stpattn1c0[i] = -999999;
01900       ntprec->trk.stpt0[i] = -999999;
01901       ntprec->trk.stpt1[i] = -999999;
01902       ntprec->trk.stptcal0t0[i] = -999999;
01903       ntprec->trk.stptcal1t0[i] = -999999;
01904     }
01905 
01906   }
01907   ntprec->trk.nstrip = nstrip;
01908 
01909 }

Generated on Sat Nov 21 22:48:00 2009 for loon by  doxygen 1.3.9.1