00001
00002
00003
00004
00005
00006
00008 #include <cmath>
00009 #include <iostream>
00010 using namespace std;
00011
00012 #include "TClonesArray.h"
00013 #include "CandNtupleSR/Module/NtpSRModule.h"
00014 #include "CandNtupleSR/Module/NtpSRBleachFiller.h"
00015 #include "CandNtupleSR/NtpSRRecord.h"
00016 #include "CandNtupleSR/NtpSRShieldStrip.h"
00017 #include "CandNtupleSR/NtpSRStrip.h"
00018 #include "CandNtupleSR/NtpSRSlice.h"
00019 #include "CandNtupleSR/NtpSRCluster.h"
00020 #include "CandNtupleSR/NtpSRShower.h"
00021 #include "CandNtupleSR/NtpSRSubShowerSummary.h"
00022 #include "CandNtupleSR/NtpSRTrack.h"
00023 #include "CandNtupleSR/NtpSREvent.h"
00024 #include "CandNtupleSR/NtpSRBleach.h"
00025 #include "CandNtupleSR/NtpSRWindow.h"
00026 #include "CandNtupleSR/NtpSRShieldExpected.h"
00027 #include "DatabaseInterface/DbiResultPtr.h"
00028 #include "MCNtuple/NtpMCRecord.h"
00029 #include "StandardNtuple/NtpStRecord.h"
00030 #include "MessageService/MsgService.h"
00031 #include "JobControl/JobCModuleRegistry.h"
00032 #include "JobControl/JobCommand.h"
00033 #include "MinosObjectMap/MomNavigator.h"
00034 #include "RawData/RawDigit.h"
00035 #include "RawData/RawRecord.h"
00036 #include "RawData/RawDaqSnarlHeader.h"
00037 #include "RawData/RawDigitDataBlock.h"
00038 #include "Record/RecCandHeader.h"
00039 #include "CandData/CandRecord.h"
00040 #include "CandData/CandHeader.h"
00041 #include "RecoBase/CandTrackListHandle.h"
00042 #include "RecoBase/CandTrackHandle.h"
00043 #include "RecoBase/CandStripListHandle.h"
00044 #include "RecoBase/CandStripHandle.h"
00045 #include "RecoBase/CandClusterListHandle.h"
00046 #include "RecoBase/CandClusterHandle.h"
00047 #include "RecoBase/CandSliceListHandle.h"
00048 #include "RecoBase/CandSliceHandle.h"
00049 #include "RecoBase/CandShowerListHandle.h"
00050 #include "RecoBase/CandShowerHandle.h"
00051 #include "RecoBase/CandEventListHandle.h"
00052 #include "RecoBase/CandEventHandle.h"
00053 #include "RecoBase/LinearFit.h"
00054 #include "RecoBase/PropagationVelocity.h"
00055 #include "CandShowerSR/CandShowerSRHandle.h"
00056 #include "CandSubShowerSR/CandSubShowerSRHandle.h"
00057 #include "CandSubShowerSR/CandSubShowerSRListHandle.h"
00058 #include "CandTrackSR/CandTrackSRHandle.h"
00059 #include "CandTrackSR/CandTrackSRListHandle.h"
00060 #include "CandFitTrackSR/CandFitTrackSRHandle.h"
00061 #include "CandFitTrackCam/CandFitTrackCamHandle.h"
00062 #include "Calibrator/Calibrator.h"
00063 #include "CandDigit/CandDeMuxDigitListHandle.h"
00064 #include "UgliGeometry/UgliGeomHandle.h"
00065 #include "Plex/PlexHandle.h"
00066 #include "Plex/PlexVetoShieldHack.h"
00067 #include "Plex/PlexPlaneId.h"
00068 #include "AstroUtil/Ast.h"
00069 #include "AstroUtil/AstTime.h"
00070 #include "AstroUtil/AstCoordinate.h"
00071 #include "Record/SimSnarlRecord.h"
00072 #include "Record/SimSnarlHeader.h"
00073 #include "CandShield/ShieldGeom.h"
00074 #include "CandShield/CandShieldSR.h"
00075 #include "CandShield/ShieldProj.h"
00076 #include "DcsUser/CoilStatus.h"
00077 #include "DcsUser/CoilTools.h"
00078 #include "DcsUser/BfldDbiCoilState.h"
00079 #include "DcsUser/DbuHvFromSingles.h"
00080 #include "DcsUser/Dcs_Mag_Near.h"
00081 #include "CandNtupleSR/NtpSRDataQuality.h"
00082 #include "CandNtupleSR/NtpSRDeadChip.h"
00083 #include "CandMorgue/CandDataQualityHandle.h"
00084 #include "CandMorgue/CandDeadChipHandle.h"
00085 #include "DataUtil/PlaneOutline.h"
00086 #include <cassert>
00087
00088 #define FIRST_ND_SPECTROMETER_PLANE 121
00089 #define FIRST_PINST_STRIP_U 0
00090 #define LAST_PINST_STRIP_U 63
00091 #define FIRST_FULL_PINST_STRIP_U 26
00092 #define LAST_FULL_PINST_STRIP_U 89
00093 #define FIRST_PINST_STRIP_V 4
00094 #define LAST_PINST_STRIP_V 67
00095
00096 ClassImp(NtpSRModule)
00097
00098
00099
00100
00101 CVSID("$Id: NtpSRModule.cxx,v 1.125 2007/10/15 03:06:47 schubert Exp $");
00102 JOBMODULE(NtpSRModule, "NtpSRModule",
00103 "A module for filling SR ntuple records.");
00104
00105 const Double_t NtpSRModule::kCos45 = 0.70710678;
00106
00107
00108
00109
00110
00111
00112
00113 const Registry& NtpSRModule::DefaultConfig() const {
00114
00115
00116
00117
00118
00119
00120
00121
00122 MSG("NtpSR",Msg::kDebug) <<
00123 "NtpSRModule::DefaultConfig" << endl;
00124
00125 static Registry r;
00126 std::string name = this->JobCModule::GetName();
00127 name += ".config.default";
00128 r.SetName(name.c_str());
00129
00130 r.UnLockValues();
00131 r.Set("PreTrigger",50.);
00132 r.Set("PostTrigger",150.);
00133 r.Set("CandRecordName", "PrimaryCandidateRecord");
00134 r.Set("SimSnarlRecordName","");
00135 r.Set("RecordName", "Primary");
00136 r.Set("RecordTitle", "Created by NtpSRModule");
00137 r.Set("HvSinglesTask",1);
00138 r.Set("UseStandard",0);
00139 r.Set("WindowPlaneExtn",10);
00140 r.Set("WindowTimeExtn",100e-9);
00141 r.LockValues();
00142
00143 return r;
00144 }
00145
00146
00147
00148 void NtpSRModule::Config(const Registry& r) {
00149
00150
00151
00152
00153
00154
00155
00156
00157 MSG("NtpSR",Msg::kDebug) << "NtpSRModule::Config" << endl;
00158
00159 Int_t tmpi;
00160 Double_t tmpd;
00161 const Char_t* tmps;
00162
00163 if ( r.Get("PreTrigger",tmpi) ) fPreTrigger = fabs((Double_t)(tmpi)*1.e-9);
00164 if ( r.Get("PreTrigger",tmpd) ) fPreTrigger = fabs(tmpd*1.e-9);
00165 if ( r.Get("PostTrigger",tmpi)) fPostTrigger = fabs((Double_t)(tmpi)*1.e-9);
00166 if ( r.Get("PostTrigger",tmpd)) fPostTrigger = fabs(tmpd*1.e-9);
00167
00168 if ( r.Get("CandRecordName",tmps) ) fCandRecordName = tmps;
00169 if ( r.Get("SimSnarlRecordName",tmps) ) fSimSnarlRecordName = tmps;
00170 if ( r.Get("RecordName", tmps) ) fRecordName = tmps;
00171 if ( r.Get("RecordTitle", tmps) ) fRecordTitle = tmps;
00172
00173 if ( r.Get("HvSinglesTask", tmpi) ) fHvSinglesTask = tmpi;
00174
00175 if ( r.Get("UseStandard",tmpi) ) fUseStandard = tmpi;
00176
00177 if ( r.Get("WindowPlaneExtn",tmpi) ) fWindowPlaneExtn = tmpi;
00178 if ( r.Get("WindowTimeExtn",tmpd) ) fWindowTimeExtn = tmpd;
00179
00180 }
00181
00182
00183
00184 JobCResult NtpSRModule::Reco(MomNavigator *mom) {
00185
00186
00187
00188
00189
00190
00191
00192
00193 JobCResult result(JobCResult::kPassed);
00194 MSG("NtpSR",Msg::kDebug) << "NtpSRModule::Reco" << endl;
00195
00196
00197 fStripUidMap.clear();
00198 fSliceUidMap.clear();
00199 fClusterUidMap.clear();
00200 fShowerUidMap.clear();
00201 fTrackUidMap.clear();
00202 fEventUidMap.clear();
00203 fUnassocStripUidMap.clear();
00204
00205
00206 assert(mom);
00207
00208
00209 const CandRecord* cndrec = dynamic_cast<const CandRecord*>
00210 (mom->GetFragment("CandRecord",fCandRecordName.c_str()));
00211 const SimSnarlRecord* simrec = dynamic_cast<const SimSnarlRecord*>
00212 (mom->GetFragment("SimSnarlRecord",fSimSnarlRecordName.c_str()));
00213 const RawRecord* rawrec = dynamic_cast<const RawRecord*>
00214 (mom->GetFragment("RawRecord","","DaqSnarl"));
00215
00216 if (!cndrec && !simrec ) {
00217 MSG("NtpSR",Msg::kWarning)
00218 << "No CandRecord of name " << fRecordName
00219 << " or SimSnarlRecord in Mom" << endl;
00220 result.SetWarning().SetFailed();
00221 return result;
00222 }
00223
00224 TClonesArray* ntpstparrptr = 0;
00225 TClonesArray* ntpexparrptr = 0;
00226 TClonesArray* ntpslcarrptr = 0;
00227 TClonesArray* ntpcluarrptr = 0;
00228 TClonesArray* ntpshwarrptr = 0;
00229 TClonesArray* ntptrkarrptr = 0;
00230 TClonesArray* ntpevtarrptr = 0;
00231 TClonesArray* ntpvetostparrptr = 0;
00232 TClonesArray* ntpdeadchipsptr = 0;
00233 NtpSREventSummary* ntpeventsummaryptr = 0;
00234 NtpSRShieldSummary* ntpshieldsummaryptr = 0;
00235 NtpSRCosmicRay* ntpcosmicrayptr = 0;
00236 NtpSRDmxStatus* ntpdmxstatusptr = 0;
00237 NtpSRDetStatus* ntpdetstatusptr = 0;
00238 NtpSRDataQuality* ntpdataqualityptr = 0;
00239
00240 NtpSRRecord* ntpsrrec = 0;
00241 NtpStRecord* ntpstrec = 0;
00242 const VldContext* ntpvldcptr = 0;
00243 if ( fUseStandard ) {
00244
00245 ntpstrec = dynamic_cast<NtpStRecord*>(mom->GetFragment("NtpStRecord",
00246 fRecordName.c_str()));
00247 if ( !ntpstrec ) {
00248 MSG("NtpSR",Msg::kWarning)
00249 << "No NtpStRecord in Mom of name " << fRecordName.c_str()
00250 << "\nMust call NtpStModule::Get first." << endl;
00251 result.SetWarning().SetFailed();
00252 return result;
00253 }
00254 if ( !cndrec ) return result;
00255 ntpstparrptr = ntpstrec -> stp;
00256 ntpslcarrptr = ntpstrec -> slc;
00257 ntpcluarrptr = ntpstrec -> clu;
00258 ntpshwarrptr = ntpstrec -> shw;
00259 ntptrkarrptr = ntpstrec -> trk;
00260 ntpevtarrptr = ntpstrec -> evt;
00261 ntpvetostparrptr = ntpstrec->vetostp;
00262 ntpexparrptr = ntpstrec->vetoexp;
00263 ntpdeadchipsptr = ntpstrec->deadchips;
00264 ntpeventsummaryptr = &(ntpstrec->evthdr);
00265 ntpshieldsummaryptr = &(ntpstrec->vetohdr);
00266 ntpcosmicrayptr = &(ntpstrec->crhdr);
00267 ntpdmxstatusptr = &(ntpstrec->dmxstatus);
00268 ntpdetstatusptr = &(ntpstrec->detstatus);
00269 ntpdataqualityptr = &(ntpstrec->dataquality);
00270 ntpvldcptr = ntpstrec->GetVldContext();
00271 }
00272 else {
00273 if ( !cndrec ) {
00274
00275 const SimSnarlHeader* simhdr = simrec->GetSimSnarlHeader();
00276
00277 RecCandHeader ntphdr(simhdr->GetVldContext(),simhdr->GetRun(),
00278 simhdr->GetSubRun(),simhdr->GetRunType(),simhdr->GetErrorCode(),
00279 simhdr->GetSnarl(),simhdr->GetTrigSrc(),
00280 simhdr->GetTimeFrame(),simhdr->GetRemoteSpillType(),-1);
00281 ntpsrrec = new NtpSRRecord(ntphdr);
00282 ntpsrrec -> SetName(fRecordName.c_str());
00283 ntpsrrec -> SetTitle(fRecordTitle.c_str());
00284 mom -> AdoptFragment(ntpsrrec);
00285 return result;
00286 }
00287
00288
00289
00290 const CandHeader* cndhdr = cndrec -> GetCandHeader();
00291 if (!rawrec) {
00292 MSG("NtpSR",Msg::kWarning) << "No DaqSnarl RawRecord in Mom"
00293 <<"\nShield data will not be filled and header will be incomplete."
00294 << endl;
00295 result.SetWarning();
00296 }
00297 if (rawrec) {
00298 const RawDaqSnarlHeader* rawhdr = dynamic_cast<const RawDaqSnarlHeader*>
00299 (rawrec -> GetRawHeader());
00300 RecCandHeader ntphdr(rawhdr->GetVldContext(),rawhdr->GetRun(),
00301 rawhdr->GetSubRun(),rawhdr->GetRunType(),rawhdr->GetErrorCode(),
00302 rawhdr->GetSnarl(),rawhdr->GetTrigSrc(),rawhdr->GetTimeFrameNum(),
00303 rawhdr->GetRemoteSpillType(),cndhdr->GetEvent());
00304 ntpsrrec = new NtpSRRecord(ntphdr);
00305 }
00306 else {
00307
00308
00309 RecCandHeader ntphdr(cndhdr->GetVldContext(),cndhdr->GetRun(),
00310 -1,-1,0,cndhdr->GetSnarl(),0,-1,-1,cndhdr->GetEvent());
00311 ntpsrrec = new NtpSRRecord(ntphdr);
00312 }
00313 ntpsrrec -> SetName(fRecordName.c_str());
00314 ntpsrrec -> SetTitle(fRecordTitle.c_str());
00315
00316 ntpstparrptr = ntpsrrec -> stp;
00317 ntpslcarrptr = ntpsrrec -> slc;
00318 ntpcluarrptr = ntpsrrec -> clu;
00319 ntpshwarrptr = ntpsrrec -> shw;
00320 ntptrkarrptr = ntpsrrec -> trk;
00321 ntpevtarrptr = ntpsrrec -> evt;
00322 ntpdeadchipsptr = ntpsrrec->deadchips;
00323 ntpvetostparrptr = ntpsrrec->vetostp;
00324 ntpexparrptr = ntpsrrec->vetoexp;
00325 ntpeventsummaryptr = &(ntpsrrec->evthdr);
00326 ntpshieldsummaryptr = &(ntpsrrec->vetohdr);
00327 ntpcosmicrayptr = &(ntpsrrec->crhdr);
00328 ntpdmxstatusptr = &(ntpsrrec->dmxstatus);
00329 ntpdetstatusptr = &(ntpsrrec->detstatus);
00330 ntpdataqualityptr = &(ntpsrrec->dataquality);
00331 ntpvldcptr = ntpsrrec->GetVldContext();
00332 }
00333
00334 TClonesArray& ntpstriparray = *(ntpstparrptr);
00335 TClonesArray& ntpslicearray = *(ntpslcarrptr);
00336 TClonesArray& ntpshieldexpected = *(ntpexparrptr);
00337 TClonesArray& ntpclusterarray = *(ntpcluarrptr);
00338 TClonesArray& ntpshowerarray = *(ntpshwarrptr);
00339 TClonesArray& ntptrackarray = *(ntptrkarrptr);
00340 TClonesArray& ntpeventarray = *(ntpevtarrptr);
00341 TClonesArray& ntpdeadchips = *(ntpdeadchipsptr);
00342 NtpSREventSummary& ntpeventsummary = *(ntpeventsummaryptr);
00343 NtpSRShieldSummary& ntpshieldsummary = *(ntpshieldsummaryptr);
00344 NtpSRCosmicRay& ntpcosmicray = *(ntpcosmicrayptr);
00345 NtpSRDmxStatus& ntpdmxstatus = *(ntpdmxstatusptr);
00346 NtpSRDetStatus& ntpdetstatus = *(ntpdetstatusptr);
00347 NtpSRDataQuality& ntpdataquality = *(ntpdataqualityptr);
00348 const VldContext& vldc = *(ntpvldcptr);
00349
00350
00351 Calibrator::Instance().Reset(vldc);
00352
00353 this -> FillNtpStrip(ntpstriparray,cndrec);
00354 this -> FillNtpSlice(ntpslicearray,cndrec);
00355 this -> FillNtpCluster(ntpclusterarray,cndrec);
00356 this -> FillNtpShower(ntpshowerarray,cndrec);
00357 this -> FillNtpTrack(ntptrackarray,cndrec);
00358 this -> FillNtpEvent(ntpeventarray,cndrec,rawrec);
00359
00360 this -> FillNtpDmxStatus(ntpdmxstatus,cndrec);
00361 this -> FillNtpDetStatus(ntpdetstatus,vldc);
00362 this -> FillNtpDataQuality(ntpdataquality,ntpdeadchips,cndrec);
00363 this -> FillNtpEventSummary(ntpeventsummary,cndrec,rawrec);
00364
00365 const NtpSRTrack* ntptrack = 0;
00366 if(ntpeventsummary.ntrack > 0) ntptrack
00367 = dynamic_cast<NtpSRTrack*>(ntptrackarray.At(0));
00368 if ( ntptrack ) this->FillNtpTrackCosmicRay(ntpcosmicray,ntptrack,vldc);
00369
00370 Detector::Detector_t dettype = vldc.GetDetector();
00371 if ( rawrec && dettype == Detector::kFar ) {
00372 TClonesArray& ntpshieldstriparray = *(ntpvetostparrptr);
00373 this->FillNtpShield(ntpshieldstriparray,ntpshieldexpected,ntpshieldsummary,ntptrack,rawrec);
00374 }
00375
00376
00377 if ( !fUseStandard ) mom -> AdoptFragment(ntpsrrec);
00378
00379 return result;
00380
00381 }
00382
00383 void NtpSRModule::FillNtpStrip(TClonesArray& ntpstriparray,
00384 const CandRecord* cndrec) {
00385
00386
00387
00388
00389
00390
00391
00392
00393 MSG("NtpSR",Msg::kDebug) << "NtpSRModule::FillNtpStrip" << endl;
00394
00395 const CandStripListHandle *striplisthandle
00396 = dynamic_cast <const CandStripListHandle*>
00397 (cndrec->FindCandHandle("CandStripListHandle"));
00398 if ( !striplisthandle ) return;
00399
00400 const CandEventListHandle *eventlisthandle
00401 = dynamic_cast <const CandEventListHandle*>
00402 (cndrec -> FindCandHandle("CandEventListHandle"));
00403
00404 Int_t nstrip = 0;
00405 TIter stripItr(striplisthandle->GetDaughterIterator());
00406 while ( CandStripHandle* strip=dynamic_cast<CandStripHandle*>(stripItr())) {
00407 Int_t uid = strip->GetUidInt();
00408 fStripUidMap.insert(std::make_pair(uid,nstrip));
00409
00410
00411 if (eventlisthandle) {
00412 Bool_t gotMatch = false;
00413 TIter eventItr(eventlisthandle->GetDaughterIterator());
00414 while(CandEventHandle* event =
00415 dynamic_cast<CandEventHandle*> (eventItr())) {
00416 TIter eventstripItr(event->GetDaughterIterator());
00417 while(CandStripHandle *eventstrip =
00418 dynamic_cast<CandStripHandle*>(eventstripItr())) {
00419 if(eventstrip->IsEquivalent(strip)) {
00420 gotMatch = true;
00421 break;
00422 }
00423 }
00424 if(gotMatch) break;
00425 }
00426 fUnassocStripUidMap.insert(std::make_pair(uid,!gotMatch));
00427 }
00428
00429
00430 NtpSRStrip* ntpstrip = new((ntpstriparray)[nstrip++])NtpSRStrip();
00431
00432 ntpstrip->index = nstrip-1;
00433 ntpstrip->planeview = strip->GetPlaneView();
00434 ntpstrip->ndigit = strip->GetNDigit();
00435 ntpstrip->demuxveto = strip->GetDemuxVetoFlag();
00436 ntpstrip->strip = strip->GetStrip();
00437 ntpstrip->plane = strip->GetPlane();
00438 ntpstrip->tpos = strip->GetTPos();
00439 ntpstrip->z = strip->GetZPos();
00440
00441
00442 bool negEndDone = false;
00443 bool posEndDone = false;
00444 TIter digitItr(strip -> GetDaughterIterator());
00445 while (CandDigitHandle* digit=dynamic_cast<CandDigitHandle*>(digitItr())) {
00446 const RawChannelId& rawch = digit->GetChannelId();
00447 Int_t pmtindex = rawch.GetCrate()*108 + rawch.GetVarcId()*36
00448 +rawch.GetVmm()*6 + rawch.GetVaAdcSel()*3+rawch.GetVaChip();
00449 if ( !negEndDone &&
00450 digit -> GetPlexSEIdAltL().GetEnd()==StripEnd::kNegative) {
00451 ntpstrip->SetPmtIndex(pmtindex,0); negEndDone = true;
00452 }
00453 else if( !posEndDone &&
00454 digit -> GetPlexSEIdAltL().GetEnd()==StripEnd::kPositive) {
00455 ntpstrip->SetPmtIndex(pmtindex,1); posEndDone = true;
00456 }
00457 }
00458
00459
00460 for (UInt_t i = 0; i < 2; i++ ) {
00461 StripEnd::EStripEnd stripend = StripEnd::kNegative;
00462 if (i == 1) stripend = StripEnd::kPositive;
00463 if ( strip->GetNDigit(stripend) > 0 ) {
00464 NtpSRPulseHeight ph;
00465 ph.raw = strip->GetCharge(CalDigitType::kNone,stripend);
00466 ph.siglin = strip->GetCharge(CalDigitType::kSigLin,stripend);
00467 ph.sigcor = strip->GetCharge(CalDigitType::kSigCorr,stripend);
00468 ph.pe = strip->GetCharge(CalDigitType::kPE,stripend);
00469 ntpstrip->SetPh(ph,i);
00470 ntpstrip->SetTime(strip->GetTime(stripend),i);
00471 }
00472 }
00473 MSG("NtpSR",Msg::kVerbose) << (*ntpstrip) << endl;
00474 }
00475
00476
00477 return;
00478
00479 }
00480
00481 void NtpSRModule::FillNtpShieldStrip(TClonesArray& ntpshieldstriparray,
00482 TClonesArray& ntpshieldexpected,
00483 NtpSRShieldSummary& ntpshieldsummary,
00484 const NtpSRTrack* ntptrack,
00485 const RawRecord* rawrec) {
00486
00487
00488
00489
00490
00491
00492
00493
00494
00495 const VldContext& vldc = *(rawrec->GetVldContext());
00496 PlexHandle plexhandle(vldc);
00497 UgliGeomHandle ugh(vldc);
00498
00499 Int_t minend = TMath::Min(StripEnd::kNegative,StripEnd::kPositive);
00500
00501 Bool_t wasExpected=false;
00502 Int_t hitstrip[2] = { 0, 0 };
00503 Double_t zproj = 0;
00504
00505 Double_t propagation_velocity = PropagationVelocity::Velocity();
00506
00507 Double_t tracktime = 0;
00508 if ( ntptrack ) tracktime = ntptrack->vtx.t;
00509
00510 PlexVetoShieldHack plexvetoshieldhack;
00511 TIter rdbit = rawrec->GetRawBlockIter();
00512 TObject* tobject;
00513
00514 std::vector<NtpSRShieldStrip*> shieldlist[3];
00515 while ((tobject = rdbit())) {
00516 RawDigitDataBlock *rdb = dynamic_cast<RawDigitDataBlock*>(tobject);
00517 if ( !rdb ) continue;
00518 TIter rdit = rdb->GetDatumIter();
00519 RawDigit* rd;
00520 while ((rd = dynamic_cast<RawDigit*>(rdit()))) {
00521 RawChannelId rawch = rd->GetChannel();
00522 PlexSEIdAltL plexaltl = plexhandle.GetSEIdAltL(rawch);
00523 Int_t stripend = 0;
00524 if ( plexaltl.IsVetoShield() ) {
00525 PlexStripEndId oldseid = plexaltl[0].GetSEId();
00526 const PlexStripEndId& newseid
00527 = plexvetoshieldhack.RenumberMuxToMdl(vldc,oldseid);
00528 stripend = 1 - (Int_t)(newseid.GetEnd()-minend);
00529 assert(stripend >= 0 && stripend < 2);
00530
00531 NtpSRShieldStrip* ntpshieldstrip = new NtpSRShieldStrip();
00532 ntpshieldstrip->ndigit++;
00533 ntpshieldstrip->pln = fShGeom->GetAssociatedPlank(newseid.GetPlane(),newseid.GetStrip(),0);
00534 ntpshieldstrip->plank = fShGeom->GetAssociatedPlank(newseid.GetPlane(),newseid.GetStrip(),1);
00535
00536
00537 wasExpected=false;
00538 NtpSRShieldExpected *ntpshexp;
00539 for(int ii=1;ii<=ntpshieldsummary.exphits;ii++){
00540 ntpshexp = dynamic_cast<NtpSRShieldExpected *>(ntpshieldexpected[ii-1]);
00541 if(ntpshexp->plane==ntpshieldstrip->pln && ntpshexp->plank==ntpshieldstrip->plank){
00542 ntpshexp->isfound=1;
00543 wasExpected=true;
00544 hitstrip[0]=ntpshexp->stripinplank[0];
00545 hitstrip[1]=ntpshexp->stripinplank[1];
00546 zproj=ntpshexp->projz;
00547 }
00548 }
00549
00550 ntpshieldstrip->adc[stripend] = rd->GetADC();
00551 ntpshieldstrip->pmtindex[stripend] = rawch.GetCrate()*108
00552 + rawch.GetVarcId()*36 + rawch.GetVmm()*6 + rawch.GetVaAdcSel()*3
00553 + rawch.GetVaChip();
00554 Int_t vach2pixel[18]={0,0,15,1,16,2,11,5,12,6,7,9,8,10,3,14,4,13};
00555 if ( rawch.GetVaChannel() >= 2 && rawch.GetVaChannel() <= 17 ) {
00556 ntpshieldstrip->pmtpixel[stripend]
00557 = vach2pixel[rawch.GetVaChannel()];
00558 }
00559 ntpshieldstrip->timeraw[stripend]
00560 = rd->GetTDC()*1.5625e-9 + 26.e-9;
00561
00562
00563 ntpshieldstrip->time[stripend] = Calibrator::Instance().GetCalibratedTime(ntpshieldstrip->timeraw[stripend],ntpshieldstrip->adc[stripend],newseid);
00564
00565 UgliStripHandle ush(ugh.GetStripHandle(newseid));
00566 if ( ush.IsValid() ) {
00567 TVector3 stripxyz0(ush.GlobalPos(-ush.GetHalfLength()));
00568 TVector3 stripxyz1(ush.GlobalPos(ush.GetHalfLength()));
00569 ntpshieldstrip -> x = fShGeom->GetPlank_X(ntpshieldstrip->pln,ntpshieldstrip->plank);
00570 ntpshieldstrip -> y = fShGeom->GetPlank_Y(ntpshieldstrip->pln,ntpshieldstrip->plank);
00571 ntpshieldstrip -> z[0] = min(stripxyz0[2],stripxyz1[2]);
00572 ntpshieldstrip -> z[1] = max(stripxyz0[2],stripxyz1[2]);
00573
00574
00575 if(wasExpected==false){
00576
00577 for ( int i = 0; i < plexaltl.GetSize(); i++ ) {
00578 PlexStripEndId stpoldseid=plexaltl[i].GetSEId();
00579 const PlexStripEndId& stpnewseid
00580 = plexvetoshieldhack.RenumberMuxToMdl(vldc,stpoldseid);
00581 UgliStripHandle stpush(ugh.GetStripHandle(stpnewseid));
00582 Int_t stpstripend
00583 = 1 - (Int_t)(stpnewseid.GetEnd()-minend);
00584 assert(stpstripend >= 0 && stpstripend < 2);
00585 ntpshieldstrip -> wlspigtail[1]
00586 += stpush.WlsPigtail(StripEnd::kNegative);
00587 ntpshieldstrip -> wlspigtail[0]
00588 += stpush.WlsPigtail(StripEnd::kPositive);
00589 ntpshieldstrip -> clearlen[1]
00590 += stpush.ClearFiber(StripEnd::kNegative);
00591 ntpshieldstrip -> clearlen[0]
00592 += stpush.ClearFiber(StripEnd::kPositive);
00593 }
00594 if ( plexaltl.GetSize() > 0 ) {
00595 for ( int iend = 0; iend < 2; iend++ ) {
00596 ntpshieldstrip->wlspigtail[iend] /=(Float_t)(plexaltl.GetSize());
00597 ntpshieldstrip->clearlen[iend] /= (Float_t)(plexaltl.GetSize());
00598 }
00599 }
00600 }else{
00601 ntpshieldstrip->wlspigtail[0]=fShGeom->GetStripWls(hitstrip[0],hitstrip[1],0);
00602 ntpshieldstrip->wlspigtail[1]=fShGeom->GetStripWls(hitstrip[0],hitstrip[1],1);
00603 ntpshieldstrip->clearlen[0]=fShGeom->GetPlaneClearFiber(hitstrip[0],0);
00604 ntpshieldstrip->clearlen[1]=fShGeom->GetPlaneClearFiber(hitstrip[0],1);
00605 }
00606
00607
00608
00609 Double_t fiberlen = ntpshieldstrip->wlspigtail[stripend]
00610 + ntpshieldstrip->clearlen[stripend];
00611 if ( fiberlen > 20. ) fiberlen = 20.;
00612 ntpshieldstrip->time[stripend] -= fiberlen/propagation_velocity;
00613 }
00614
00615
00616 if(wasExpected==false){
00617 if ( ntptrack && ntpshieldsummary.projz >= ntpshieldstrip->z[0] &&
00618 ntpshieldsummary.projz <= ntpshieldstrip->z[1] ) {
00619 ntpshieldstrip->time[stripend] -= fabs(ntpshieldsummary.projz
00620 -ntpshieldstrip->z[stripend])/propagation_velocity;
00621 }
00622 } else{
00623 if(stripend==0){
00624 ntpshieldstrip->time[stripend] -= fabs(zproj-fShGeom->GetPlank_Z(ntpshieldstrip->pln,ntpshieldstrip->plank)+4.)/propagation_velocity;
00625 }else if(stripend==1){
00626 ntpshieldstrip->time[stripend] -= fabs(4.+fShGeom->GetPlank_Z(ntpshieldstrip->pln,ntpshieldstrip->plank)-zproj)/propagation_velocity;
00627 }
00628 }
00629
00630
00631 if( ntptrack ){
00632 if((fShGeom->WhatSection(ntpshieldstrip->pln)==fShGeom->ClosestTwoSections(ntptrack->vtx.z,0) && 1.e9*fabs(ntptrack->vtx.t-ntpshieldstrip->time[stripend])<100) || (fShGeom->WhatSection(ntpshieldstrip->pln)==fShGeom->ClosestTwoSections(ntptrack->vtx.z,1) && 1.e9*fabs(ntptrack->vtx.t-ntpshieldstrip->time[stripend])<100)){
00633 ntpshieldsummary.found2sect=true;
00634 }
00635 }
00636
00637
00638
00639 Int_t ishieldtime = 1;
00640 if ( (ntpshieldstrip->time[stripend]-tracktime)< -1.*fabs(fPreTrigger))
00641 ishieldtime = 0;
00642 else if((ntpshieldstrip->time[stripend]-tracktime)>fabs(fPostTrigger))
00643 ishieldtime = 2;
00644 ntpshieldsummary.ndigit[ishieldtime]++;
00645 ntpshieldsummary.adc[ishieldtime] += rd->GetADC();;
00646
00647
00648
00649
00650 bool isFound = false;
00651 std::vector<NtpSRShieldStrip*>::iterator vsItr;
00652 for ( vsItr = shieldlist[ishieldtime].begin();
00653 vsItr!= shieldlist[ishieldtime].end() && !isFound; vsItr++) {
00654 NtpSRShieldStrip* liststrip = *vsItr;
00655 if ( liststrip->pln == ntpshieldstrip->pln &&
00656 liststrip->plank == ntpshieldstrip->plank ) {
00657 isFound = true;
00658 liststrip->ndigit += ntpshieldstrip->ndigit;
00659 if ( liststrip->adc[stripend] == 0 ) {
00660 liststrip->timeraw[stripend] = ntpshieldstrip->timeraw[stripend];
00661 liststrip->time[stripend] = ntpshieldstrip->time[stripend];
00662 liststrip->wlspigtail[stripend]
00663 = ntpshieldstrip->wlspigtail[stripend];
00664 liststrip->clearlen[stripend]=ntpshieldstrip->clearlen[stripend];
00665 liststrip->pmtindex[stripend]=ntpshieldstrip->pmtindex[stripend];
00666 liststrip->pmtpixel[stripend]=ntpshieldstrip->pmtpixel[stripend];
00667 }
00668 liststrip->adc[stripend] += ntpshieldstrip->adc[stripend];
00669 delete ntpshieldstrip; ntpshieldstrip = 0;
00670 }
00671 }
00672 if ( !isFound ) {
00673 shieldlist[ishieldtime].push_back(ntpshieldstrip);
00674 ntpshieldsummary.nplank[ishieldtime]++;
00675 }
00676 }
00677 }
00678
00679
00680
00681
00682
00683 Int_t nshieldstrip = 0;
00684 for ( int ishieldtime = 0; ishieldtime < 3; ishieldtime++ ) {
00685 std::vector<NtpSRShieldStrip*>::iterator vsItr;
00686 for ( vsItr = shieldlist[ishieldtime].begin();
00687 vsItr!= shieldlist[ishieldtime].end(); vsItr++) {
00688 NtpSRShieldStrip* liststrip = *vsItr;
00689
00690 NtpSRShieldStrip& ntpstrip =
00691 *(new(ntpshieldstriparray[nshieldstrip++])NtpSRShieldStrip(*liststrip));
00692 ntpstrip.index = nshieldstrip - 1;
00693
00694
00695 NtpSRShieldExpected *ntpshexp;
00696 for(int ii=1;ii<=ntpshieldsummary.exphits;ii++){
00697 ntpshexp = dynamic_cast<NtpSRShieldExpected *>(ntpshieldexpected[ii-1]);
00698 if(ntpshexp->plane==ntpstrip.pln && ntpshexp->plank==ntpstrip.plank){
00699 ntpshexp->stripdigit=ntpstrip.index;
00700 }
00701 }
00702
00703 if ( ntptrack ) {
00704
00705
00706 const NtpSRVertex& vtx = ntptrack -> vtx;
00707
00708 Float_t dx;
00709 Int_t sign=1;
00710 ShieldProj sp(vtx.x,vtx.y,vtx.z,vtx.dcosx,vtx.dcosy,vtx.dcosz,ntpstrip.pln,ntpstrip.plank,fShGeom);
00711 if(sp.GetProjInter_Z() >= ntpstrip.z[0] && sp.GetProjInter_Z() <= ntpstrip.z[1]){
00712 if(fShGeom->IsVertical(ntpstrip.pln)){
00713 if(sp.GetProjInter_Y()-(fShGeom->GetPlank_Y(ntpstrip.pln,ntpstrip.plank)) < 0){
00714 sign=-1;
00715 }
00716 } else{
00717 if(sp.GetProjInter_X()-(fShGeom->GetPlank_X(ntpstrip.pln,ntpstrip.plank)) < 0){
00718 sign=-1;
00719 }
00720 }
00721 dx=sp.GetProjDis()*sign;
00722 for(Int_t ie = 0; ie<2; ie++){
00723 if(fabs(dx)<fabs(ntpshieldsummary.dx[ishieldtime])){
00724 ntpshieldsummary.dx[ishieldtime] = dx;
00725 ntpshieldsummary.dxvetostp[ishieldtime] = nshieldstrip-1;
00726 }
00727 }
00728 }
00729 }
00730 delete liststrip; liststrip = 0;
00731 MSG("NtpSR",Msg::kDebug) << ntpstrip << endl;
00732 }
00733 }
00734 }
00735
00736 return;
00737
00738 }
00739
00740 void NtpSRModule::FillNtpTrackProjectionToShield(NtpSRShieldSummary&
00741 ntpshieldsummary, const NtpSRTrack* ntptrack, const CandShieldSR& cssh) {
00742
00743
00744
00745
00746
00747
00748
00749
00750
00751
00752
00753
00754
00755 MSG("NtpSR",Msg::kVerbose) << "NtpSRModule::FillNtpTrackProjectionToShield"
00756 << endl;
00757
00758 if ( ! ntptrack ) return;
00759
00760 const NtpSRVertex& vtx = ntptrack->vtx;
00761
00762 ntpshieldsummary.ishit = 0;
00763 if(cssh.IsVetoHit() == true){
00764 ntpshieldsummary.ishit = 1;
00765 }
00766 ntpshieldsummary.projx=cssh.GetCandShieldInter_X();
00767 ntpshieldsummary.projy=cssh.GetCandShieldInter_Y();
00768 ntpshieldsummary.projz=cssh.GetCandShieldInter_Z();
00769
00770 ntpshieldsummary.exphits=cssh.HitsInShield();
00771
00772 Float_t xzproj[2] = {-999.,-999.};
00773 Int_t ishit=0;
00774 if ( vtx.dcosy != 0. ) {
00775 xzproj[0] = vtx.x + vtx.dcosx/vtx.dcosy*(4.37-vtx.y);
00776 xzproj[1] = vtx.z + vtx.dcosz/vtx.dcosy*(4.37-vtx.y);
00777 if (fabs(xzproj[0])<2.03 && (xzproj[1]-vtx.z)*vtx.dcosz<0.) {
00778 if (xzproj[1]>0.0&&xzproj[1]<30.) ishit = 1;
00779 if ((xzproj[1]>0.0&&xzproj[1]<30.) || !ishit) {
00780 ntpshieldsummary.dcos = vtx.dcosy;
00781 }
00782 }
00783 xzproj[0] = vtx.x + vtx.dcosx/vtx.dcosy*(3.06-vtx.y);
00784 xzproj[1] = vtx.z + vtx.dcosz/vtx.dcosy*(3.06-vtx.y);
00785 if (fabs(xzproj[0])>3.34 && (xzproj[1]-vtx.z)*vtx.dcosz<0.) {
00786 if (fabs(xzproj[0])<6.05&&xzproj[1]>0.0&&xzproj[1]<30.)
00787 ishit = 1;
00788 if ((xzproj[1]>0.0&&xzproj[1]<30.) || !ishit) {
00789 ntpshieldsummary.dcos = vtx.dcosy;
00790 }
00791 }
00792 }
00793 if (vtx.dcosu!=0.) {
00794 xzproj[0] = vtx.v + vtx.dcosv/vtx.dcosu*(4.45-vtx.u);
00795 xzproj[1] = vtx.z + vtx.dcosz/vtx.dcosu*(4.45-vtx.u);
00796 if (xzproj[0]>-.22624&&xzproj[0]<1.5554
00797 && (xzproj[1]-vtx.z)*vtx.dcosz<0.) {
00798 if (xzproj[1]>0.0&&xzproj[1]<30.) ishit = 1;
00799 if ((xzproj[1]>0.0&&xzproj[1]<30.) || !ishit) {
00800 ntpshieldsummary.dcos = vtx.dcosu;
00801 }
00802 }
00803 }
00804 if (vtx.dcosv!=0.) {
00805 xzproj[0] = vtx.u +vtx.dcosu/vtx.dcosv*(4.45-vtx.v);
00806 xzproj[1] = vtx.z + vtx.dcosz/vtx.dcosv*(4.45-vtx.v);
00807 if (xzproj[0]>-.22624&&xzproj[0]<1.5554
00808 && (xzproj[1]-vtx.z)*vtx.dcosz<0.) {
00809 if (xzproj[1]>0.0&&xzproj[1]<30.) ishit = 1;
00810 if ((xzproj[1]>0.0&&xzproj[1]<30.) || !ishit) {
00811 ntpshieldsummary.dcos = vtx.dcosv;
00812 }
00813 }
00814 }
00815 if (vtx.dcosx!=0.) {
00816 xzproj[0] = vtx.y + vtx.dcosy/vtx.dcosx*(6.7-vtx.x);
00817 xzproj[1] = vtx.z + vtx.dcosz/vtx.dcosx*(6.7-vtx.x);
00818 if (xzproj[0]>1.20&&xzproj[0]<4.64 && (xzproj[1]-vtx.z)*vtx.dcosz<0.) {
00819 if (xzproj[1]>0.0&&xzproj[1]<30.) ishit = 1;
00820 if ((xzproj[1]>0.0&&xzproj[1]<30.) || !ishit) {
00821 ntpshieldsummary.dcos = vtx.dcosx;
00822 }
00823 }
00824 xzproj[0] = vtx.y + vtx.dcosy/vtx.dcosx*(-6.7-vtx.x);
00825 xzproj[1] = vtx.z + vtx.dcosz/vtx.dcosx*(-6.7-vtx.x);
00826 if (xzproj[0]>1.20&&xzproj[0]<4.64 && (xzproj[1]-vtx.z)*vtx.dcosz<0.) {
00827 if (xzproj[1]>0.0&&xzproj[1]<30.) ishit = 1;
00828 if ((xzproj[1]>0.0&&xzproj[1]<30.) || !ishit) {
00829 ntpshieldsummary.dcos = vtx.dcosx;
00830 }
00831 }
00832 xzproj[0] = vtx.y + vtx.dcosy/vtx.dcosx*(4.3-vtx.x);
00833 xzproj[1] = vtx.z + vtx.dcosz/vtx.dcosx*(4.3-vtx.x);
00834 if (xzproj[0]>-1.5&&xzproj[0]<0.5 && (xzproj[1]-vtx.z)*vtx.dcosz<0.) {
00835 if (xzproj[1]>0.0&&xzproj[1]<30.) ishit = 1;
00836 if ((xzproj[1]>0.0&&xzproj[1]<30.) || !ishit) {
00837 ntpshieldsummary.dcos = vtx.dcosx;
00838 }
00839 }
00840 xzproj[0] = vtx.y + vtx.dcosy/vtx.dcosx*(-4.3-vtx.x);
00841 xzproj[1] = vtx.z + vtx.dcosz/vtx.dcosx*(-4.3-vtx.x);
00842 if (xzproj[0]>-1.5&&xzproj[0]<0.5 && (xzproj[1]-vtx.z)*vtx.dcosz<0.) {
00843 if (xzproj[1]>0.0&&xzproj[1]<30.) ishit = 1;
00844 if ((xzproj[1]>0.0&&xzproj[1]<30.) || !ishit) {
00845 ntpshieldsummary.dcos = vtx.dcosx;
00846 }
00847 }
00848 }
00849
00850 return;
00851
00852 }
00853
00854 void NtpSRModule::FillNtpShieldExpected(TClonesArray& ntpshieldexpected,
00855 const CandShieldSR& cssh) {
00856
00857
00858
00859
00860
00861
00862
00863
00864
00865
00866 MSG("NtpSR",Msg::kVerbose) << "NtpSRModule::FillNtpShieldExpected"
00867 << endl;
00868
00869 Int_t cont=0;
00870 for(int ii=1;ii<=cssh.HitsInShield();ii++){
00871 cont+=1;
00872 NtpSRShieldExpected* ntpshexp = new((ntpshieldexpected)[cont-1])NtpSRShieldExpected();
00873 ntpshexp->plane=cssh.GetCandShieldPlane(ii);
00874 ntpshexp->plank=cssh.GetCandShieldStrip0(ii);
00875 ntpshexp->stripinplank[0]=cssh.GetStripInPlank(ii,0);
00876 ntpshexp->stripinplank[1]=cssh.GetStripInPlank(ii,1);
00877 ntpshexp->projx=cssh.GetCandShieldInter_X(ii);
00878 ntpshexp->projy=cssh.GetCandShieldInter_Y(ii);
00879 ntpshexp->projz=cssh.GetCandShieldInter_Z(ii);
00880 ntpshexp->centerdis=cssh.GetInterCenterDis(ii);
00881 ntpshexp->index=cont-1;
00882 ntpshexp->isfound=0;
00883 ntpshexp->stripdigit=-1;
00884
00885 }
00886 }
00887
00888 void NtpSRModule::FillNtpShield(TClonesArray& ntpshieldstriparray,
00889 TClonesArray& ntpshieldexpected,
00890 NtpSRShieldSummary& ntpshieldsummary,
00891 const NtpSRTrack* ntptrack,
00892 const RawRecord* rawrec) {
00893
00894
00895
00896
00897
00898
00899
00900
00901
00902
00903 MSG("NtpSR",Msg::kDebug) << "NtpSRModule::FillNtpShield" << endl;
00904
00905
00906 const VldContext& vldc = *(rawrec->GetVldContext());
00907 if( !fShGeom ){
00908 fShGeom = new ShieldGeom(vldc);
00909 }
00910 else{
00911 fShGeom->Reinitialize(vldc);
00912 }
00913 if( ntptrack ){
00914 const NtpSRVertex& vtx = ntptrack->vtx;
00915 const NtpSRVertex& end = ntptrack->end;
00916 const CandShieldSR cssh(vtx,end,fShGeom);
00917
00918
00919
00920
00921 this -> FillNtpTrackProjectionToShield(ntpshieldsummary,ntptrack,cssh);
00922 this -> FillNtpShieldExpected(ntpshieldexpected,cssh);
00923 }
00924 this -> FillNtpShieldStrip(ntpshieldstriparray,ntpshieldexpected,ntpshieldsummary,ntptrack,rawrec);
00925
00926 MSG("NtpSR",Msg::kDebug) << ntpshieldsummary << endl;
00927
00928 return;
00929
00930 }
00931
00932 void NtpSRModule::FillNtpSlice(TClonesArray& ntpslicearray,
00933 const CandRecord* cndrec) {
00934
00935
00936
00937
00938
00939
00940
00941
00942 MSG("NtpSR",Msg::kDebug) << "NtpSRModule::FillNtpSlice" << endl;
00943
00944 const CandSliceListHandle *slicelisthandle
00945 = dynamic_cast <const CandSliceListHandle*>
00946 (cndrec -> FindCandHandle("CandSliceListHandle"));
00947 if ( !slicelisthandle ) return;
00948
00949 Int_t nslice = 0;
00950 TIter sliceItr(slicelisthandle->GetDaughterIterator());
00951 while ( CandSliceHandle* slice=dynamic_cast<CandSliceHandle*>(sliceItr())) {
00952 Int_t uid = slice->GetUidInt();
00953
00954 fSliceUidMap.insert(std::make_pair(uid,nslice));
00955
00956 NtpSRSlice* ntpslice = new(ntpslicearray[nslice++])
00957 NtpSRSlice(slice->GetNStrip());
00958 ntpslice->index = nslice-1;
00959 ntpslice->ndigit = slice->GetNDigit();
00960
00961 TIter slicestripItr(slice->GetDaughterIterator());
00962 Int_t nslicestrip = 0;
00963 while ( CandStripHandle *slicestrip = dynamic_cast<CandStripHandle*>
00964 (slicestripItr())) {
00965 Int_t uid = slicestrip->GetUidInt();
00966 std::map<int,int>::iterator uidItr;
00967 uidItr = fStripUidMap.find(uid);
00968 if ( uidItr == fStripUidMap.end() ) {
00969 MSG("NtpSR",Msg::kError)
00970 << "Slice strip w/Uid " << uid
00971 << " does not match any in strip list."
00972 << "\n slc stp entry will not be properly filled." << endl;
00973 }
00974 else {
00975 Int_t stripindex = uidItr->second;
00976 ntpslice->AddStripAt(stripindex,nslicestrip);
00977 }
00978 nslicestrip++;
00979 }
00980
00981
00982 ntpslice->ph.raw = slice->GetCharge(CalDigitType::kNone);
00983 ntpslice->ph.siglin = slice->GetCharge(CalDigitType::kSigLin);
00984 ntpslice->ph.sigcor = slice->GetCharge(CalDigitType::kSigCorr);
00985 ntpslice->ph.pe = slice->GetCharge(CalDigitType::kPE);
00986
00987
00988 ntpslice->plane.n =slice->GetNPlane();
00989 ntpslice->plane.beg = slice->GetBegPlane();
00990 ntpslice->plane.end = slice->GetEndPlane();
00991 ntpslice->plane.nu = slice->GetNPlane(PlaneView::kU);
00992 ntpslice->plane.begu = slice->GetBegPlane(PlaneView::kU);
00993 ntpslice->plane.endu = slice->GetEndPlane(PlaneView::kU);
00994 ntpslice->plane.nv = slice->GetNPlane(PlaneView::kV);
00995 ntpslice->plane.begv = slice->GetBegPlane(PlaneView::kV);
00996 ntpslice->plane.endv = slice->GetEndPlane(PlaneView::kV);
00997
00998 MSG("NtpSR",Msg::kDebug) << "CandSlice uid " << slice-> GetUidInt()
00999 << "\n" << (*ntpslice) << endl;
01000 }
01001
01002 return;
01003 }
01004
01005 void NtpSRModule::FillNtpCluster(TClonesArray& ntpclusterarray,
01006 const CandRecord* cndrec) {
01007
01008
01009
01010
01011
01012
01013
01014
01015 MSG("NtpSR",Msg::kDebug) << "NtpSRModule::FillNtpCluster" << endl;
01016
01017 Int_t ncluster = 0;
01018
01019 const CandSubShowerSRListHandle *clusterlisthandle
01020 = dynamic_cast <const CandSubShowerSRListHandle*>
01021 (cndrec -> FindCandHandle("CandSubShowerSRListHandle"));
01022
01023 if ( !clusterlisthandle ) {
01024
01025
01026 const CandClusterListHandle *clulisthandle
01027 = dynamic_cast <const CandClusterListHandle*>
01028 (cndrec -> FindCandHandle("CandClusterListHandle"));
01029 if( !clulisthandle ) return;
01030
01031 TIter cluItr(clulisthandle->GetDaughterIterator());
01032 while (CandClusterHandle* clus=dynamic_cast<CandClusterHandle*>(cluItr())){
01033 Int_t uid = clus->GetUidInt();
01034 fClusterUidMap.insert(std::make_pair(uid,ncluster));
01035
01036
01037 NtpSRCluster* ntpcluster
01038 = new(ntpclusterarray[ncluster++])NtpSRCluster(clus->GetNStrip());
01039
01040
01041 ntpcluster->index = ncluster - 1;
01042
01043
01044 const CandSliceHandle* clusterslice = clus -> GetCandSlice();
01045 if ( clusterslice ) {
01046 std::map<int,int>::iterator uidItr
01047 = fSliceUidMap.find(clusterslice->GetUidInt());
01048 if ( uidItr == fSliceUidMap.end() ) {
01049 MSG("NtpSR",Msg::kError)
01050 << "Cluster slice w/Uid " << clusterslice->GetUidInt()
01051 << " does not match any in slice list."
01052 << "\n clu slc will not be properly filled." << endl;
01053 }
01054 else {
01055 ntpcluster->slc = uidItr->second;
01056 }
01057 }
01058 else {
01059 MSG("NtpSR",Msg::kWarning) << "No associated Slice found for cluster"
01060 << endl;
01061 ntpcluster->slc = -1;
01062 }
01063
01064 TIter clusterstripItr(clus->GetDaughterIterator());
01065 Int_t nclusterstrip = 0;
01066 while ( CandStripHandle *clusterstrip = dynamic_cast<CandStripHandle*>
01067 (clusterstripItr()) ) {
01068 std::map<int,int>::iterator uidItr
01069 = fStripUidMap.find(clusterstrip->GetUidInt());
01070 if ( uidItr == fStripUidMap.end() ) {
01071 MSG("NtpSR",Msg::kError)
01072 << "Cluster strip w/Uid " << clusterstrip->GetUidInt()
01073 << " does not match any in strip list."
01074 << "\n clu stp entry will not be properly filled." << endl;
01075 }
01076 else {
01077 Int_t stripindex = uidItr->second;
01078 ntpcluster->AddStripAt(stripindex,nclusterstrip);
01079 }
01080 nclusterstrip++;
01081 }
01082
01083 ntpcluster->planeview = clus->GetPlaneView();
01084 ntpcluster->nplane = clus->GetNPlane();
01085 ntpcluster->begplane = clus->GetBegPlane();
01086 ntpcluster->endplane = clus->GetEndPlane();
01087 ntpcluster->id = ClusterType::kUnknown;
01088
01089
01090 ntpcluster->ph.raw = clus->GetCharge();
01091 ntpcluster->ph.siglin = 0.;
01092 ntpcluster->ph.sigcor = 0.;
01093 ntpcluster->ph.pe = 0.;
01094 ntpcluster->ph.sigmap = 0.;
01095 ntpcluster->ph.mip = 0.;
01096 ntpcluster->ph.gev = 0.;
01097 MSG("NtpSR",Msg::kDebug) << "CandCluster uid "
01098 << clus -> GetUidInt() << "\n"
01099 << (*ntpcluster) << endl;
01100 }
01101 return;
01102 }
01103
01104
01105 CandSubShowerSRHandleItr clusterItr(clusterlisthandle->GetDaughterIterator());
01106 CandSubShowerSRHandleKeyFunc *engKF = clusterItr.CreateKeyFunc();
01107 engKF->SetFun(CandSubShowerSRHandle::KeyFromViewEnergy);
01108 clusterItr.GetSet()->AdoptSortKeyFunc(engKF);
01109 engKF = 0;
01110 while (CandSubShowerSRHandle* cluster=dynamic_cast<CandSubShowerSRHandle*>
01111 (clusterItr())) {
01112 Int_t uid = cluster->GetUidInt();
01113 fClusterUidMap.insert(std::make_pair(uid,ncluster));
01114
01115
01116 NtpSRCluster* ntpcluster
01117 = new(ntpclusterarray[ncluster++])NtpSRCluster(cluster->GetNStrip());
01118
01119
01120 ntpcluster->index = ncluster - 1;
01121
01122
01123 const CandSliceHandle* clusterslice = cluster -> GetCandSlice();
01124 if ( clusterslice ) {
01125 std::map<int,int>::iterator uidItr
01126 = fSliceUidMap.find(clusterslice->GetUidInt());
01127 if ( uidItr == fSliceUidMap.end() ) {
01128 MSG("NtpSR",Msg::kError)
01129 << "Cluster slice w/Uid " << clusterslice->GetUidInt()
01130 << " does not match any in slice list."
01131 << "\n clu slc will not be properly filled." << endl;
01132 }
01133 else {
01134 ntpcluster->slc = uidItr->second;
01135 }
01136 }
01137 else {
01138 MSG("NtpSR",Msg::kWarning) << "No associated Slice found for cluster"
01139 << endl;
01140 ntpcluster->slc = -1;
01141 }
01142
01143 ntpcluster->ndigit = cluster->GetNDigit();
01144
01145 TIter clusterstripItr(cluster->GetDaughterIterator());
01146 Int_t nclusterstrip = 0;
01147 while ( CandStripHandle *clusterstrip = dynamic_cast<CandStripHandle*>
01148 (clusterstripItr()) ) {
01149 std::map<int,int>::iterator uidItr
01150 = fStripUidMap.find(clusterstrip->GetUidInt());
01151 if ( uidItr == fStripUidMap.end() ) {
01152 MSG("NtpSR",Msg::kError)
01153 << "Cluster strip w/Uid " << clusterstrip->GetUidInt()
01154 << " does not match any in strip list."
01155 << "\n clu stp entry will not be properly filled." << endl;
01156 }
01157 else {
01158 Int_t stripindex = uidItr->second;
01159 ntpcluster->AddStripAt(stripindex,nclusterstrip);
01160 }
01161 nclusterstrip++;
01162 }
01163
01164 ntpcluster->planeview = cluster->GetPlaneView();
01165 ntpcluster->nplane = cluster->GetNPlane();
01166 ntpcluster->begplane = cluster->GetBegPlane();
01167 ntpcluster->endplane = cluster->GetEndPlane();
01168
01169 ntpcluster->zvtx = cluster->GetVtxZ();
01170 if(ntpcluster->planeview==PlaneView::kU ||
01171 ntpcluster->planeview==PlaneView::kX) ntpcluster->tposvtx = cluster->GetVtxU();
01172 else if(ntpcluster->planeview==PlaneView::kV ||
01173 ntpcluster->planeview==PlaneView::kY) ntpcluster->tposvtx = cluster->GetVtxV();
01174
01175 ntpcluster->slope = cluster->GetSlope();
01176 ntpcluster->avgdev = cluster->GetAvgDev();
01177 ntpcluster->id = cluster->GetClusterID();
01178 ntpcluster->probem = cluster->GetProbEM();
01179
01180
01181 ntpcluster->ph.raw = cluster->GetCharge(CalStripType::kNone);
01182 ntpcluster->ph.siglin = cluster->GetCharge(CalStripType::kSigLin);
01183 ntpcluster->ph.sigcor = cluster->GetCharge(CalStripType::kSigCorr);
01184 ntpcluster->ph.pe = cluster->GetCharge(CalStripType::kPE);
01185 ntpcluster->ph.sigmap = cluster->GetCharge(CalStripType::kSigMapped);
01186 ntpcluster->ph.mip = cluster->GetCharge(CalStripType::kMIP);
01187 ntpcluster->ph.gev = cluster->GetEnergy();
01188 MSG("NtpSR",Msg::kDebug) << "CandSubShowerSR uid "
01189 << cluster -> GetUidInt() << "\n"
01190 << (*ntpcluster) << endl;
01191 }
01192
01193 return;
01194 }
01195
01196 void NtpSRModule::FillNtpShower(TClonesArray& ntpshowerarray,
01197 const CandRecord* cndrec) {
01198
01199
01200
01201
01202
01203
01204
01205
01206 MSG("NtpSR",Msg::kDebug) << "NtpSRModule::FillNtpShower" << endl;
01207
01208 const CandShowerListHandle *showerlisthandle
01209 = dynamic_cast <const CandShowerListHandle*>
01210 (cndrec -> FindCandHandle("CandShowerListHandle"));
01211 if ( !showerlisthandle ) return;
01212
01213 const CandTrackListHandle *tracklisthandle
01214 = dynamic_cast <const CandTrackListHandle*>
01215 (cndrec -> FindCandHandle("CandTrackListHandle"));
01216
01217 Int_t nshower = 0;
01218 TIter showerItr(showerlisthandle->GetDaughterIterator());
01219 while (CandShowerHandle* shower=dynamic_cast<CandShowerHandle*>
01220 (showerItr())) {
01221 Int_t uid = shower->GetUidInt();
01222 fShowerUidMap.insert(std::make_pair(uid,nshower));
01223
01224 NtpSRShower* ntpshower = 0;
01225 if (shower -> InheritsFrom("CandShowerSRHandle")) {
01226 CandShowerSRHandle* showerSR = dynamic_cast<CandShowerSRHandle*>(shower);
01227 ntpshower
01228 = new(ntpshowerarray[nshower++])NtpSRShower(showerSR->GetNStrip(),
01229 showerSR->GetNumSubShowersU()+showerSR->GetNumSubShowersV());
01230 }
01231 else {
01232 ntpshower=new(ntpshowerarray[nshower++])NtpSRShower(shower->GetNStrip(),
01233 shower->GetLastCluster()+1);
01234 }
01235 ntpshower -> nstpcnt = shower->GetNStrip();
01236
01237
01238 ntpshower->index = nshower - 1;
01239
01240
01241 const CandSliceHandle* showerslice = shower -> GetCandSlice();
01242 if ( showerslice ) {
01243 std::map<int,int>::iterator uidItr;
01244 uidItr = fSliceUidMap.find(showerslice->GetUidInt());
01245 if ( uidItr == fSliceUidMap.end() ) {
01246 MSG("NtpSR",Msg::kError)
01247 << "Shower slice w/Uid "
01248 << showerslice->GetUidInt()
01249 << " does not match any in slice list."
01250 << "\n shw slc will not be properly filled." << endl;
01251 }
01252 else {
01253 ntpshower->slc = uidItr->second;
01254 }
01255 }
01256 else {
01257 MSG("NtpSR",Msg::kWarning) << "No associated Slice found for shower"
01258 << endl;
01259 ntpshower->slc = -1;
01260 }
01261
01262 ntpshower->ndigit = shower->GetNDigit();
01263 TIter showerstripItr(shower->GetDaughterIterator());
01264 Int_t nshowerstrip = 0;
01265 while ( CandStripHandle *showerstrip = dynamic_cast<CandStripHandle*>
01266 (showerstripItr()) ) {
01267 Int_t uid = showerstrip->GetUidInt();
01268 std::map<int,int>::iterator uidItr;
01269 uidItr = fStripUidMap.find(uid);
01270 if ( uidItr == fStripUidMap.end() ) {
01271 MSG("NtpSR",Msg::kError)
01272 << "Shower strip w/Uid " << uid
01273 << " does not match any in strip list."
01274 << "\n shw stp entry will not be properly filled." << endl;
01275 }
01276 else {
01277 Int_t stripindex = uidItr->second;
01278 ntpshower->AddStripAt(stripindex,nshowerstrip);
01279 }
01280
01281
01282 Int_t iplane = showerstrip->GetPlane();
01283 ntpshower->stpu[nshowerstrip] = shower->GetU(iplane);
01284 ntpshower->stpv[nshowerstrip] = shower->GetV(iplane);
01285 ntpshower->stpx[nshowerstrip] = kCos45*(ntpshower->stpu[nshowerstrip]
01286 -ntpshower->stpv[nshowerstrip]);
01287 ntpshower->stpy[nshowerstrip] = kCos45*(ntpshower->stpu[nshowerstrip]
01288 +ntpshower->stpv[nshowerstrip]);
01289 ntpshower->stpz[nshowerstrip] = showerstrip->GetZPos();
01290
01291 for (UInt_t iend = 0; iend < 2; iend++ ) {
01292 StripEnd::EStripEnd stripend = StripEnd::kNegative;
01293 if (iend == 1) stripend = StripEnd::kPositive;
01294 if ( showerstrip->GetNDigit(stripend) > 0 ) {
01295 Float_t sigmap = shower->GetStripCharge(showerstrip,
01296 CalStripType::kSigMapped,stripend);
01297 Float_t mip = shower->GetStripCharge(showerstrip,
01298 CalStripType::kMIP,stripend);
01299 Float_t gev = shower->GetStripCharge(showerstrip,
01300 CalStripType::kGeV,stripend);
01301 ntpshower->SetPh(sigmap,mip,gev,nshowerstrip,iend);
01302 ntpshower->SetTime(shower->GetT(iplane,stripend),nshowerstrip,iend);
01303 double c0 = Calibrator::Instance().DecalStripToStrip(1.0,
01304 showerstrip->GetStripEndId(stripend));
01305 ntpshower->SetAttnC0(c0,nshowerstrip,iend);
01306
01307 double t0 = Calibrator::Instance().DecalTime(0.0,
01308 showerstrip->GetCharge(CalDigitType::kNone,stripend),
01309 showerstrip->GetStripEndId(stripend));
01310 ntpshower->SetCalT0(t0,nshowerstrip,iend);
01311 }
01312 }
01313 nshowerstrip++;
01314 }
01315
01316
01317 ntpshower->ncluster = 0;
01318 ntpshower->nUcluster = 0;
01319 ntpshower->nVcluster = 0;
01320
01321 if (shower->InheritsFrom("CandShowerSRHandle")) {
01322 CandShowerSRHandle* showerSR = dynamic_cast<CandShowerSRHandle*>(shower);
01323 ntpshower->ncluster = showerSR->GetLastSubShower()+1;
01324 if ( showerSR->GetLastSubShower()+1 > 0 ) {
01325 ntpshower->nUcluster = showerSR->GetNumSubShowersU();
01326 ntpshower->nVcluster = showerSR->GetNumSubShowersV();
01327 for (int i = 0; i < ntpshower->ncluster; i++) {
01328 const CandSubShowerSRHandle* showercluster=showerSR->GetSubShower(i);
01329 Int_t uid = showercluster->GetUidInt();
01330 std::map<int,int>::iterator uidItr = fClusterUidMap.find(uid);
01331 if ( uidItr == fClusterUidMap.end() ) {
01332 MSG("NtpSR",Msg::kError)
01333 << "Shower CandSubShowerSR w/Uid " << uid
01334 << " does not match any in cluster list."
01335 << " shw clu entry will not be properly filled." << endl;
01336 }
01337 else {
01338 Int_t clusterindex = uidItr->second;
01339 ntpshower -> AddClusterAt(clusterindex,i);
01340 }
01341 }
01342 this->FillNtpSubShowerSummary(ntpshower,showerSR,tracklisthandle);
01343 }
01344 }
01345 else {
01346 ntpshower->ncluster = shower->GetLastCluster()+1;
01347 for (int i = 0; i < ntpshower->ncluster; i++ ) {
01348 const CandClusterHandle* showercluster = shower->GetCluster(i);
01349 Int_t uid = showercluster->GetUidInt();
01350 std::map<int,int>::iterator uidItr = fClusterUidMap.find(uid);
01351 if ( uidItr == fClusterUidMap.end() ) {
01352 MSG("NtpSR",Msg::kError)
01353 << "\nShower CandCluster of index " << i << " (of "
01354 << ntpshower->ncluster << " clusters) w/Uid " << uid
01355 << " does not match any in cluster list."
01356 << "\nshw clu entry for CandShower of index "
01357 << ntpshower->index << " w/Uid " << shower->GetUidInt()
01358 << " will not be properly filled!" << endl;
01359 }
01360 else {
01361 Int_t clusterindex = uidItr->second;
01362 ntpshower -> AddClusterAt(clusterindex,i);
01363 }
01364 if ( showercluster->GetPlaneView() == PlaneView::kU )
01365 ntpshower->nUcluster += 1;
01366 else if (showercluster->GetPlaneView() == PlaneView::kV )
01367 ntpshower->nVcluster += 1;
01368 }
01369 }
01370
01371
01372 ntpshower->plane.n = shower->GetNPlane();
01373 ntpshower->plane.nu = shower->GetNPlane(PlaneView::kU);
01374 ntpshower->plane.nv = shower->GetNPlane(PlaneView::kV);
01375 ntpshower->plane.beg = shower->GetBegPlane();
01376 ntpshower->plane.begu = shower->GetBegPlane(PlaneView::kU);
01377 ntpshower->plane.begv = shower->GetBegPlane(PlaneView::kV);
01378 ntpshower->plane.end = shower->GetEndPlane();
01379 ntpshower->plane.endu = shower->GetEndPlane(PlaneView::kU);
01380 ntpshower->plane.endv = shower->GetEndPlane(PlaneView::kV);
01381 ntpshower->contained = shower->IsContained();
01382
01383 ntpshower->ph.raw = shower->GetCharge(CalStripType::kNone);
01384 ntpshower->ph.siglin = shower->GetCharge(CalStripType::kSigLin);
01385 ntpshower->ph.sigcor = shower->GetCharge(CalStripType::kSigCorr);
01386 ntpshower->ph.pe = shower->GetCharge(CalStripType::kPE);
01387 ntpshower->ph.sigmap = shower->GetCharge(CalStripType::kSigMapped);
01388 ntpshower->ph.mip = shower->GetCharge(CalStripType::kMIP);
01389 ntpshower->ph.gev = shower->GetEnergy(CandShowerHandle::kWtCC);
01390 ntpshower->shwph.wtCCgev = shower->GetEnergy(CandShowerHandle::kWtCC);
01391 ntpshower->shwph.linCCgev = shower->GetEnergy(CandShowerHandle::kCC);
01392 ntpshower->shwph.wtNCgev = shower->GetEnergy(CandShowerHandle::kWtNC);
01393 ntpshower->shwph.linNCgev = shower->GetEnergy(CandShowerHandle::kNC);
01394 ntpshower->shwph.EMgev = shower->GetEnergy(CandShowerHandle::kEM);
01395
01396
01397 NtpSRVertex& vtx = ntpshower->vtx;
01398 vtx.u = shower->GetVtxU();
01399 vtx.v = shower->GetVtxV();
01400 vtx.x = kCos45*(vtx.u-vtx.v);
01401 vtx.y = kCos45*(vtx.u+vtx.v);
01402 vtx.z = shower->GetVtxZ();
01403 vtx.t = shower->GetVtxT();
01404 vtx.plane = shower->GetVtxPlane();
01405 vtx.dcosu = shower->GetDirCosU();
01406 vtx.dcosv = shower->GetDirCosV();
01407 vtx.dcosx = kCos45*(vtx.dcosu-vtx.dcosv);
01408 vtx.dcosy = kCos45*(vtx.dcosu+vtx.dcosv);
01409 vtx.dcosz = shower->GetDirCosZ();
01410 MSG("NtpSR",Msg::kDebug) << "CandShower uid "
01411 << shower -> GetUidInt() << "\n"
01412 << (*ntpshower) << endl;
01413 }
01414
01415 return;
01416 }
01417
01418
01419 void NtpSRModule::FillNtpEvent(TClonesArray& ntpeventarray,
01420 const CandRecord* cndrec,
01421 const RawRecord* rawrec) {
01422
01423
01424
01425
01426
01427
01428
01429
01430 MSG("NtpSR",Msg::kDebug) << "NtpSRModule::FillNtpEvent" << endl;
01431
01432 const CandEventListHandle *eventlisthandle
01433 = dynamic_cast <const CandEventListHandle*>
01434 (cndrec -> FindCandHandle("CandEventListHandle"));
01435 if (!eventlisthandle) return;
01436
01437 Int_t nevent = 0;
01438 TIter eventItr(eventlisthandle->GetDaughterIterator());
01439 while (CandEventHandle* event=dynamic_cast<CandEventHandle*>
01440 (eventItr())) {
01441 Int_t uid = event->GetUidInt();
01442 fEventUidMap.insert(std::make_pair(uid,nevent));
01443
01444 NtpSREvent* ntpevent
01445 = new(ntpeventarray[nevent++])NtpSREvent(event->GetNStrip(),
01446 event->GetLastShower()+1,event->GetLastTrack()+1);
01447 ntpevent->index = nevent-1;
01448
01449 const CandSliceHandle* eventslice = event -> GetCandSlice();
01450
01451 if ( eventslice ) {
01452 std::map<int,int>::iterator uidItr;
01453 uidItr = fSliceUidMap.find(eventslice->GetUidInt());
01454 if ( uidItr == fSliceUidMap.end() ) {
01455 MSG("NtpSR",Msg::kError)
01456 << "Event slice w/Uid "
01457 << eventslice->GetUidInt()
01458 << " does not match any in slice list."
01459 << "\n evt slc will not be properly filled." << endl;
01460 }
01461 else {
01462 ntpevent->slc = uidItr->second;
01463 }
01464 }
01465 else {
01466 MSG("NtpSR",Msg::kWarning) << "No associated Slice found for event"
01467 << endl;
01468 ntpevent->slc = -1;
01469 }
01470
01471 ntpevent->ndigit = event->GetNDigit();
01472
01473 TIter eventstripItr(event->GetDaughterIterator());
01474 Int_t neventstrip = 0;
01475
01476 std::map<int,int>::iterator uidItr;
01477 while ( CandStripHandle *eventstrip = dynamic_cast<CandStripHandle*>
01478 (eventstripItr())) {
01479 uidItr = fStripUidMap.find(eventstrip->GetUidInt());
01480 if ( uidItr == fStripUidMap.end() ) {
01481 MSG("NtpSR",Msg::kError)
01482 << "Event strip w/Uid "
01483 << eventstrip->GetUidInt()
01484 << " does not match any in strip list."
01485 << "\nevt stp array will not be properly filled." << endl;
01486 }
01487 else {
01488 Int_t stripindex = uidItr->second;
01489 ntpevent -> AddStripAt(stripindex,neventstrip);
01490 }
01491
01492 for (UInt_t iend = 0; iend < 2; iend++ ) {
01493 StripEnd::EStripEnd stripend = StripEnd::kNegative;
01494 if (iend == 1) stripend = StripEnd::kPositive;
01495 if ( eventstrip->GetNDigit(stripend) > 0 ) {
01496 Float_t sigmap = event->GetStripCharge(eventstrip,
01497 CalStripType::kSigMapped,stripend);
01498 Float_t mip = event->GetStripCharge(eventstrip,
01499 CalStripType::kMIP,stripend);
01500 Float_t gev = event->GetStripCharge(eventstrip,
01501 CalStripType::kGeV,stripend);
01502 ntpevent->SetPh(sigmap,mip,gev,neventstrip,iend);
01503 }
01504 }
01505
01506 neventstrip++;
01507 }
01508
01509 for (Int_t i = 0; i <= event->GetLastTrack(); i++ ) {
01510 const CandTrackHandle* track = event->GetTrack(i);
01511 uidItr = fTrackUidMap.find(track->GetUidInt());
01512 if ( uidItr == fTrackUidMap.end() ) {
01513 MSG("NtpSR",Msg::kError)
01514 << "Event track w/Uid "
01515 << track->GetUidInt() << " does not match any in track list."
01516 << "\nevt trk array will not be properly filled." << endl;
01517 }
01518 else {
01519 Int_t trackindex = uidItr->second;
01520 ntpevent -> AddTrackAt(trackindex,i);
01521 }
01522 }
01523 for (Int_t i = 0; i <= event->GetLastShower(); i++ ) {
01524 const CandShowerHandle* shower = event->GetShower(i);
01525 uidItr = fShowerUidMap.find(shower->GetUidInt());
01526 if ( uidItr == fShowerUidMap.end() ) {
01527 MSG("NtpSR",Msg::kError)
01528 << "Event shower w/Uid "
01529 << shower->GetUidInt() << " does not match any in shower list."
01530 << "\nevt shw array will not be properly filled." << endl;
01531 }
01532 else {
01533 Int_t showerindex = uidItr->second;
01534 ntpevent -> AddShowerAt(showerindex,i);
01535 }
01536 }
01537 ntpevent->primtrk = event->GetPrimaryTrackIndex();
01538 ntpevent->primshw = event->GetPrimaryShowerIndex();
01539
01540
01541 ntpevent->plane.n = event->GetNPlane();
01542 ntpevent->plane.nu = event->GetNPlane(PlaneView::kU);
01543 ntpevent->plane.nv = event->GetNPlane(PlaneView::kV);
01544 ntpevent->plane.beg = event->GetBegPlane();
01545 ntpevent->plane.begu = event->GetBegPlane(PlaneView::kU);
01546 ntpevent->plane.begv = event->GetBegPlane(PlaneView::kV);
01547 ntpevent->plane.end = event->GetEndPlane();
01548 ntpevent->plane.endu = event->GetEndPlane(PlaneView::kU);
01549 ntpevent->plane.endv = event->GetEndPlane(PlaneView::kV);
01550 ntpevent->contained = event->IsContained();
01551
01552 ntpevent->ph.raw = event->GetCharge(CalStripType::kNone);
01553 ntpevent->ph.siglin = event->GetCharge(CalStripType::kSigLin);
01554 ntpevent->ph.sigcor = event->GetCharge(CalStripType::kSigCorr);
01555 ntpevent->ph.pe = event->GetCharge(CalStripType::kPE);
01556 ntpevent->ph.sigmap = event->GetCharge(CalStripType::kSigMapped);
01557 ntpevent->ph.mip = event->GetCharge(CalStripType::kMIP);
01558
01559 ntpevent->ph.gev = event->GetEnergy();
01560
01561
01562 NtpSRVertex& vtx = ntpevent->vtx;
01563 vtx.u = event->GetVtxU();
01564 vtx.v = event->GetVtxV();
01565 vtx.x = kCos45*(vtx.u - vtx.v);
01566 vtx.y = kCos45*(vtx.u + vtx.v);
01567 vtx.z = event->GetVtxZ();
01568 vtx.t = event->GetVtxT();
01569 vtx.plane = event->GetVtxPlane();
01570 vtx.dcosu = event->GetVtxDirCosU();
01571 vtx.dcosv = event->GetVtxDirCosV();
01572 vtx.dcosx = kCos45*(vtx.dcosu - vtx.dcosv);
01573 vtx.dcosy = kCos45*(vtx.dcosu + vtx.dcosv);
01574 vtx.dcosz = event->GetVtxDirCosZ();
01575
01576 NtpSRVertex& end = ntpevent->end;
01577 end.u = event->GetEndU();
01578 end.v = event->GetEndV();
01579 end.x = kCos45*(end.u - end.v);
01580 end.y = kCos45*(end.u + end.v);
01581 end.z = event->GetEndZ();
01582 end.t = event->GetEndT();
01583 end.plane = event->GetEndPlane();
01584 end.dcosu = event->GetEndDirCosU();
01585 end.dcosv = event->GetEndDirCosV();
01586 end.dcosx = kCos45*(end.dcosu - end.dcosv);
01587 end.dcosy = kCos45*(end.dcosu + end.dcosv);
01588 end.dcosz = event->GetEndDirCosZ();
01589
01590 NtpSRBleach& bleach = ntpevent->bleach;
01591 FillNtpBleach(bleach,event,cndrec,rawrec);
01592
01593 FillNtpWindow(ntpevent,cndrec);
01594
01595 MSG("NtpSR",Msg::kDebug) << (*ntpevent) << endl;
01596 }
01597
01598
01599 return;
01600 }
01601
01602 void NtpSRModule::FillNtpTrack(TClonesArray& ntptrackarray,
01603 const CandRecord* cndrec) {
01604
01605
01606
01607
01608
01609
01610
01611
01612
01613 MSG("NtpSR",Msg::kDebug) << "NtpSRModule::FillNtpTrack" << endl;
01614
01615 const CandTrackListHandle *tracklisthandle
01616 = dynamic_cast <const CandTrackListHandle*>
01617 (cndrec -> FindCandHandle("CandFitTrackListHandle"));
01618 if ( !tracklisthandle ) {
01619 tracklisthandle = dynamic_cast <const CandTrackListHandle*>
01620 (cndrec -> FindCandHandle("CandTrackListHandle"));
01621 }
01622 if ( !tracklisthandle ) return;
01623
01624 CandTrackSRListHandle* tracksrlisthandle=dynamic_cast<CandTrackSRListHandle*>
01625 (cndrec->FindCandHandle("CandTrackSRListHandle"));
01626
01627 const VldContext& vld = *(cndrec->GetVldContext());
01628
01629 TIter trackItr(tracklisthandle->GetDaughterIterator());
01630 Int_t ntrack = 0;
01631 while (CandTrackHandle* track = dynamic_cast<CandTrackHandle*>(trackItr())){
01632 if(track->GetNDaughters()==0 &&
01633 track->InheritsFrom("CandFitTrackHandle") &&
01634 dynamic_cast<CandFitTrackHandle*>(track)->GetFinderTrack()){
01635 track = dynamic_cast<CandFitTrackHandle*>(track)->GetFinderTrack();
01636 }
01637 Int_t uid = track->GetUidInt();
01638 fTrackUidMap.insert(std::make_pair(uid,ntrack));
01639
01640 NtpSRTrack* ntptrack
01641 = new(ntptrackarray[ntrack++])NtpSRTrack(track->GetNStrip());
01642 ntptrack->index = ntrack-1;
01643
01644
01645 const CandSliceHandle* trackslice = track -> GetCandSlice();
01646 if ( trackslice ) {
01647 std::map<int,int>::iterator uidItr;
01648 uidItr = fSliceUidMap.find(trackslice->GetUidInt());
01649 if ( uidItr == fSliceUidMap.end() ) {
01650 MSG("NtpSR",Msg::kError)
01651 << "Track slice w/Uid "
01652 << trackslice->GetUidInt()
01653 << " does not match any in slice list."
01654 << "\n trk slc will not be properly filled." << endl;
01655 }
01656 else {
01657 ntptrack->slc = uidItr->second;
01658 }
01659 }
01660 else {
01661 MSG("NtpSR",Msg::kWarning) << "No associated Slice found for track"
01662 << endl;
01663 ntptrack->slc = -1;
01664 }
01665
01666 ntptrack->ndigit = track->GetNDigit();
01667 CandFitTrackSRHandle *fittracksr
01668 =dynamic_cast<CandFitTrackSRHandle*>(track);
01669 CandFitTrackHandle *fittrack
01670 =dynamic_cast<CandFitTrackHandle*>(track);
01671
01672
01673 NtpSRTrackPlane& plane = ntptrack->plane;
01674 plane.n = track->GetNPlane();
01675 plane.nu = track->GetNPlane(PlaneView::kU);
01676 plane.nv = track->GetNPlane(PlaneView::kV);
01677 plane.beg = track->GetBegPlane();
01678 plane.begu = track->GetBegPlane(PlaneView::kU);
01679 plane.begv = track->GetBegPlane(PlaneView::kV);
01680 plane.end = track->GetEndPlane();
01681 plane.endu = track->GetEndPlane(PlaneView::kU);
01682 plane.endv = track->GetEndPlane(PlaneView::kV);
01683 plane.ntrklike = track->GetNTrackPlane();
01684
01685
01686 NtpSRStripPulseHeight& ph = ntptrack->ph;
01687 ph.raw = track->GetCharge(CalStripType::kNone);
01688 ph.siglin = track->GetCharge(CalStripType::kSigLin);
01689 ph.sigcor = track->GetCharge(CalStripType::kSigCorr);
01690 ph.pe = track->GetCharge(CalStripType::kPE);
01691 ph.sigmap = track->GetCharge(CalStripType::kSigMapped);
01692 ph.mip = track->GetCharge(CalStripType::kMIP);
01693 ph.gev = track->GetCharge(CalStripType::kGeV);
01694
01695 ntptrack->ds = track->GetdS();
01696 ntptrack->range = track->GetRange();
01697
01698 if (tracksrlisthandle) ntptrack->cputime = tracksrlisthandle->GetCPUTime();
01699
01700 this->FillNtpTrackVertex(ntptrack,track);
01701 this->FillNtpTrackEnd(ntptrack,track);
01702 this->FillNtpTrackLinearFit(ntptrack,track);
01703
01704 this->FillNtpTrackFidVtx(ntptrack,track,vld);
01705 this->FillNtpTrackFidEnd(ntptrack,track,vld);
01706 this->FillNtpTrackFidAll(ntptrack,track,vld);
01707
01708 this->FillNtpTrackFit(ntptrack,track);
01709 this->FillNtpTrackMomentum(ntptrack,track);
01710 this->FillNtpTrackTime(ntptrack,track);
01711 ntptrack->contained = track->IsContained();
01712
01713 NtpSRCosmicRay& ntpcosmicray = ntptrack->cr;
01714 this->FillNtpTrackCosmicRay(ntpcosmicray,ntptrack,vld);
01715
01716
01717 TIter trackstripItr(track->GetDaughterIterator());
01718 Int_t ntrackstrip = 0;
01719
01720 while ( CandStripHandle *trackstrip = dynamic_cast<CandStripHandle*>
01721 (trackstripItr())) {
01722 Int_t uid = trackstrip->GetUidInt();
01723 std::map<int,int>::iterator uidItr;
01724 uidItr = fStripUidMap.find(uid);
01725 if ( uidItr == fStripUidMap.end() ) {
01726 MSG("NtpSR",Msg::kError)
01727 << "Track strip w/Uid " << uid
01728 << " does not match any in strip list."
01729 << "\n trk stp entry will not be properly filled." << endl;
01730 }
01731 else {
01732 Int_t stripindex = uidItr->second;
01733 ntptrack->AddStripAt(stripindex,ntrackstrip);
01734 }
01735
01736
01737 Int_t iplane = trackstrip->GetPlane();
01738 ntptrack->stpu[ntrackstrip] = track->GetU(iplane);
01739 ntptrack->stpv[ntrackstrip] = track->GetV(iplane);
01740 ntptrack->stpx[ntrackstrip] = kCos45*(ntptrack->stpu[ntrackstrip]
01741 -ntptrack->stpv[ntrackstrip]);
01742 ntptrack->stpy[ntrackstrip] = kCos45*(ntptrack->stpu[ntrackstrip]
01743 +ntptrack->stpv[ntrackstrip]);
01744 ntptrack->stpz[ntrackstrip] = track->GetZ(iplane);
01745
01746 ntptrack->stpds[ntrackstrip] = track->GetdS(iplane);
01747
01748
01749
01750 if ( fittrack ) {
01751 ntptrack->stpfitchi2[ntrackstrip] = fittrack->GetPlaneChi2(iplane);
01752 ntptrack->stpfitqp[ntrackstrip] = fittrack->GetPlaneQP(iplane);
01753 }
01754 else{
01755 ntptrack->stpfitchi2[ntrackstrip] = 0;
01756 ntptrack->stpfitqp[ntrackstrip] = 0;
01757 }
01758
01759 if ( fittracksr ) {
01760 ntptrack->stpfitprechi2[ntrackstrip]
01761 = fittracksr->GetPlanePreChi2(iplane);
01762 if ( fittracksr->GetBadFit(iplane) ) ntptrack->stpfit[ntrackstrip] = 0;
01763 }
01764 else{
01765 ntptrack->stpfitprechi2[ntrackstrip] =0;
01766 ntptrack->stpfit[ntrackstrip] = 0;
01767 }
01768
01769
01770 for (UInt_t i = 0; i < 2; i++ ) {
01771 StripEnd::EStripEnd stripend = StripEnd::kNegative;
01772 if (i == 1) stripend = StripEnd::kPositive;
01773 if ( trackstrip->GetNDigit(stripend) > 0 ) {
01774 Float_t sigmap = track->GetStripCharge(trackstrip,
01775 CalStripType::kSigMapped,stripend);
01776 Float_t mip = track->GetStripCharge(trackstrip,
01777 CalStripType::kMIP,stripend);
01778 Float_t gev = track->GetStripCharge(trackstrip,
01779 CalStripType::kGeV,stripend);
01780 ntptrack->SetPh(sigmap,mip,gev,ntrackstrip,i);
01781 ntptrack->SetTime(track->GetT(iplane,stripend),ntrackstrip,i);
01782
01783
01784
01785
01786
01787 double c0 = Calibrator::Instance().DecalStripToStrip(1.0,
01788 trackstrip->GetStripEndId(stripend));
01789 ntptrack->SetAttnC0(c0,ntrackstrip,i);
01790
01791 double t0 = Calibrator::Instance().DecalTime(0.0,
01792 trackstrip->GetCharge(CalDigitType::kNone,stripend),
01793 trackstrip->GetStripEndId(stripend));
01794 ntptrack->SetCalT0(t0,ntrackstrip,i);
01795 }
01796 }
01797 ntrackstrip++;
01798 }
01799
01800 MSG("NtpSR",Msg::kDebug) << (*ntptrack) << endl;
01801
01802 }
01803
01804 return;
01805
01806 }
01807
01808 void NtpSRModule::FillNtpFiducialDistance(NtpSRFiducial& fid,
01809 const NtpSRVertex& vtx,
01810 const VldContext& vld) {
01811
01812
01813
01814
01815
01816
01817
01818
01819 UgliGeomHandle ugh(vld);
01820 Float_t zextent[2];
01821 ugh.GetZExtent(zextent[0],zextent[1]);
01822 Detector::Detector_t Detector = vld.GetDetector();
01823 PlaneOutline pl;
01824 PlexPlaneId plnid(Detector,vtx.plane,false);
01825 float dist=-1; float xedge=0; float yedge=0;
01826 if(plnid.IsValid()){
01827 pl.DistanceToEdge(vtx.x, vtx.y,
01828 plnid.GetPlaneView(),
01829 plnid.GetPlaneCoverage(),
01830 dist, xedge, yedge);
01831 }
01832 fid.dr=dist;
01833 if ( fid.dr < 0. ) fid.dr = 0;
01834 fid.dz = min(vtx.z-zextent[0],zextent[1]-vtx.z);
01835 }
01836
01837 void NtpSRModule::FillNtpTrackMomentum(NtpSRTrack* ntptrack,
01838 const CandTrackHandle* track) {
01839
01840
01841
01842
01843
01844
01845
01846
01847
01848 MSG("NtpSR",Msg::kVerbose) << "NtpSRModule::FillNtpTrackMomentum" << endl;
01849 const CandFitTrackHandle *fittrack
01850 =dynamic_cast<const CandFitTrackHandle*>(track);
01851
01852 NtpSRMomentum& momentum = ntptrack->momentum;
01853 if ( fittrack ) {
01854
01855 momentum.range = fittrack->GetMomentumRange();
01856 momentum.best = fittrack->GetMomentum();
01857 Double_t pcurve = max(fittrack->GetMomentumCurve(),0.001);
01858
01859 momentum.qp = fittrack->GetEMCharge()/pcurve;
01860 momentum.qp_rangebiased = 0;
01861
01862 const CandFitTrackCamHandle *fittrackcam
01863 =dynamic_cast<const CandFitTrackCamHandle*>(track);
01864
01865 if(fittrackcam) momentum.qp_rangebiased = fittrackcam->GetRangeBiasedQP();
01866
01867 momentum.eqp = fittrack->GetVtxQPError();
01868 }
01869 else{
01870 momentum.range = track->GetMomentum();
01871 momentum.best = track->GetMomentum();
01872 }
01873
01874 return;
01875
01876 }
01877
01878 void NtpSRModule::FillNtpTrackCosmicRay(NtpSRCosmicRay& ntpcosmicray,
01879 const NtpSRTrack* ntptrack,
01880 const VldContext& vldc) {
01881
01882
01883
01884
01885
01886
01887
01888
01889
01890 MSG("NtpSR",Msg::kVerbose) << "NtpSRModule::FillNtpTrackCosmicRay" << endl;
01891 if ( !ntptrack ) return;
01892
01893 Detector::Detector_t dettype = vldc.GetDetector();
01894
01895
01896 Float_t dcosx = -(ntptrack->vtx.dcosx);
01897 Float_t dcosy = -(ntptrack->vtx.dcosy);
01898 Float_t dcosz = -(ntptrack->vtx.dcosz);
01899
01900 Double_t altitude,azimuth;
01901 AstUtil::LocalToHorizon(dcosx,dcosy,dcosz,dettype,altitude,azimuth);
01902
01903 ntpcosmicray.zenith = 90. - altitude;
01904 ntpcosmicray.azimuth = azimuth;
01905
01906 UInt_t dyear,dmonth,dday,dhour,dminute,dsec;
01907 vldc.GetTimeStamp().GetDate(kTRUE,0,&dyear,&dmonth,&dday);
01908 vldc.GetTimeStamp().GetTime(kTRUE,0,&dhour,&dminute,&dsec);
01909
01910 int year = dyear;
01911 int month = dmonth;
01912 int day = dday;
01913 double hour = dhour + (dminute + (dsec)/60.)/60.;
01914
01915 AstUtil::CalendarToJulian(year,month,day,hour,ntpcosmicray.juliandate);
01916 double longitude = AstUtil::GetDetLongitude(dettype);
01917 double latitude = AstUtil::GetDetLatitude(dettype);
01918
01919 double gast;
01920 AstUtil::JulianToGAST(ntpcosmicray.juliandate,gast);
01921 double last;
01922 AstUtil::GSTToLST(gast,longitude,last);
01923
01924 ntpcosmicray.locsiderialtime = last;
01925 double hourangle, declination;
01926 AstUtil::HorizonToEquatorial(altitude,azimuth,latitude,hourangle,
01927 declination);
01928 ntpcosmicray.dec = declination;
01929 double ra;
01930 AstUtil::EquatorialToCelestial(hourangle,gast,longitude,ra);
01931 ntpcosmicray.ra = ra;
01932
01933 ntpcosmicray.rahourangle = ntpcosmicray.ra*12./180.;
01934
01935 return;
01936
01937 }
01938
01939 void NtpSRModule::FillNtpTrackTime(NtpSRTrack* ntptrack,
01940 const CandTrackHandle* track) {
01941
01942
01943
01944
01945
01946
01947
01948
01949
01950 MSG("NtpSR",Msg::kVerbose) << "NtpSRModule::FillNtpTrackTime" << endl;
01951
01952 NtpSRTrackTime& time = ntptrack->time;
01953
01954 time.ndigit = track->GetNTimeFitDigit();
01955 time.chi2 = track->GetTimeFitChi2();
01956
01957 time.cdtds = fabs(track->GetTimeSlope())*3.e8;
01958 time.dtds = track->GetTimeSlope();
01959 time.t0 = track->GetTimeOffset();
01960
01961 time.forwardRMS = track->GetTimeForwardFitRMS();
01962 time.forwardNDOF = track->GetTimeForwardFitNDOF();
01963 time.backwardRMS = track->GetTimeBackwardFitRMS();
01964 time.backwardNDOF = track->GetTimeBackwardFitNDOF();
01965
01966 Double_t totuph[2] = {0}, totvph[2] = {0};
01967 Int_t ndutimespace[1000] = {0}, ndvtimespace[1000] = {0};
01968 Double_t dutimespace[1000] = {0}, dvtimespace[1000] = {0};
01969
01970 TIter trackstripItr(track->GetDaughterIterator());
01971 while (CandStripHandle* trackstrip = dynamic_cast<CandStripHandle*>
01972 (trackstripItr())) {
01973 Int_t iplane = trackstrip -> GetPlane();
01974 assert(iplane >= 0 && iplane < 1000);
01975 Double_t ph0 = trackstrip -> GetCharge(StripEnd::kNegative);
01976 Double_t ph1 = trackstrip -> GetCharge(StripEnd::kPositive);
01977 Double_t t0 = track->GetT(iplane,StripEnd::kNegative);
01978 Double_t t1 = track->GetT(iplane,StripEnd::kPositive);
01979
01980 if ( trackstrip -> GetPlaneView() == PlaneView::kU ) {
01981 if ( trackstrip -> GetNDigit(StripEnd::kNegative) > 0 ) {
01982 totuph[0] += ph0;
01983 time.u0 += ph0*t0;
01984 }
01985 if ( trackstrip -> GetNDigit(StripEnd::kPositive) > 0 ) {
01986 totuph[1] += ph1;
01987 time.u1 += ph1*t1;
01988 }
01989 if ( trackstrip->GetNDigit(StripEnd::kNegative) > 0 &&
01990 trackstrip->GetNDigit(StripEnd::kPositive) > 0 ) {
01991 Double_t dapos = .5*(t0-t1)*PropagationVelocity::VelocityPosErr();
01992 if (!ndutimespace[iplane] || fabs(dapos)<fabs(dutimespace[iplane])) {
01993 dutimespace[iplane] = dapos;
01994 ndutimespace[iplane] = 1;
01995 }
01996 }
01997 }
01998
01999 if ( trackstrip -> GetPlaneView() == PlaneView::kV ) {
02000 if ( trackstrip -> GetNDigit(StripEnd::kNegative) > 0 ) {
02001 totvph[0] += ph0;
02002 time.v0 += ph0*t0;
02003 }
02004 if ( trackstrip -> GetNDigit(StripEnd::kPositive) > 0 ) {
02005 totvph[1] += ph1;
02006 time.v1 += ph1*t1;
02007 }
02008 if ( trackstrip -> GetNDigit(StripEnd::kNegative) > 0 &&
02009 trackstrip -> GetNDigit(StripEnd::kPositive) > 0 ) {
02010 Double_t dapos = .5*(t0-t1)*PropagationVelocity::VelocityPosErr();
02011 if (!ndvtimespace[iplane] || fabs(dapos) < fabs(dvtimespace[iplane])) {
02012 dvtimespace[iplane] = dapos;
02013 ndvtimespace[iplane] = 1;
02014 }
02015 }
02016 }
02017 }
02018
02019 Int_t ndu = 0;
02020 Int_t ndv = 0;
02021 for (int ip = 0; ip < 1000; ip++ ) {
02022 if ( ndutimespace[ip] > 0 ) {
02023 ndu++;
02024 time.du += dutimespace[ip];
02025 }
02026 if ( ndvtimespace[ip] > 0 ) {
02027 ndv++;
02028 time.dv += dvtimespace[ip];
02029 }
02030 }
02031 if ( totuph[0] > 0. ) time.u0 /= totuph[0];
02032 if ( totvph[0] > 0. ) time.v0 /= totvph[0];
02033 if ( totuph[1] > 0. ) time.u1 /= totuph[1];
02034 if ( totvph[1] > 0. ) time.v1 /= totvph[1];
02035 if ( ndu ) time.du /= (Double_t)ndu;
02036 if ( ndv ) time.dv /= (Double_t)ndv;
02037
02038 return;
02039
02040 }
02041
02042 void NtpSRModule::FillNtpTrackFit(NtpSRTrack* ntptrack,
02043 const CandTrackHandle* track) {
02044
02045
02046
02047
02048
02049
02050
02051
02052
02053 MSG("NtpSR",Msg::kVerbose) << "NtpSRModule::FillNtpTrackFit" << endl;
02054
02055 NtpSRFitTrack& fit = ntptrack->fit;
02056
02057 const CandFitTrackHandle* fittrack
02058 = dynamic_cast<const CandFitTrackHandle*>(track);
02059 if ( fittrack ) {
02060 fit.chi2 = fittrack->GetChi2();
02061 fit.pass = fittrack->GetPass();
02062 fit.ndof = fittrack->GetNDOF();
02063 fit.niterate = fittrack->GetNIterate();
02064 fit.cputime = fittrack->GetCPUTime();
02065 fit.nswimfail = fittrack->GetNSwimFail();
02066 fit.bave = fittrack->GetBave();
02067 }
02068
02069 return;
02070
02071 }
02072
02073 void NtpSRModule::FillNtpTrackFidVtx(NtpSRTrack* ntptrack,
02074 const CandTrackHandle* track,
02075 const VldContext& vld) {
02076
02077
02078
02079
02080
02081
02082
02083
02084
02085
02086 MSG("NtpSR",Msg::kVerbose) << "NtpSRModule::FillNtpTrackFidVtx" << endl;
02087
02088 NtpSRFiducial& fid = ntptrack->fidvtx;
02089 NtpSRVertex& vtx = ntptrack->vtx;
02090 this->FillNtpFiducialDistance(fid,vtx,vld);
02091
02092 fid.trace = track->GetVtxTrace();
02093 fid.tracez = track->GetVtxTraceZ();
02094 fid.nplane = track->GetVtxnActiveUpstream();
02095 return;
02096 }
02097
02098 void NtpSRModule::FillNtpTrackFidEnd(NtpSRTrack* ntptrack,
02099 const CandTrackHandle* track,
02100 const VldContext& vld) {
02101
02102
02103
02104
02105
02106
02107
02108
02109
02110
02111 MSG("NtpSR",Msg::kVerbose) << "NtpSRModule::FillNtpTrackFidEnd" << endl;
02112
02113 NtpSRFiducial& fid = ntptrack->fidend;
02114 NtpSRVertex& end = ntptrack->end;
02115 this->FillNtpFiducialDistance(fid,end,vld);
02116
02117 fid.trace = track->GetEndTrace();
02118 fid.tracez = track->GetEndTraceZ();
02119 fid.nplane = track->GetEndnActiveDownstream();
02120 return;
02121
02122 }
02123
02124 void NtpSRModule::FillNtpTrackFidAll(NtpSRTrack* ntptrack,
02125 const CandTrackHandle* track,
02126 const VldContext& vld) {
02127
02128
02129
02130
02131
02132
02133
02134
02135
02136
02137
02138
02139 MSG("NtpSR",Msg::kVerbose) << "NtpSRModule::FillNtpFiducialAll" << endl;
02140
02141 NtpSRFiducial& fidall = ntptrack->fidall;
02142
02143 NtpSRFiducial fidtmp;
02144 NtpSRVertex vtxtmp;
02145 fidall.dr = min(ntptrack->fidvtx.dr,ntptrack->fidend.dr);
02146 fidall.dz = min(ntptrack->fidvtx.dz,ntptrack->fidend.dz);
02147 fidall.trace = min(ntptrack->fidvtx.trace,ntptrack->fidend.trace);
02148 fidall.tracez = min(ntptrack->fidvtx.tracez,ntptrack->fidend.tracez);
02149 fidall.nplane = min(ntptrack->fidvtx.nplane,ntptrack->fidend.nplane);
02150
02151 TIter trackstripItr(track->GetDaughterIterator());
02152 while(CandStripHandle* trackstrip = dynamic_cast<CandStripHandle*>
02153 (trackstripItr())) {
02154
02155 Int_t iplane = trackstrip->GetPlane();
02156 vtxtmp.plane = iplane;
02157 vtxtmp.u = track->GetU(iplane);
02158 vtxtmp.v = track->GetV(iplane);
02159 vtxtmp.x = kCos45*(vtxtmp.u-vtxtmp.v);
02160 vtxtmp.y = kCos45*(vtxtmp.u+vtxtmp.v);
02161 vtxtmp.z = track->GetZ(iplane);
02162 this->FillNtpFiducialDistance(fidtmp,vtxtmp,vld);
02163 fidall.dr = min(fidtmp.dr,fidall.dr);
02164 fidall.dz = min(fidtmp.dz,fidall.dz);
02165 }
02166
02167 return;
02168
02169 }
02170
02171 void NtpSRModule::FillNtpTrackVertex(NtpSRTrack* ntptrack,
02172 const CandTrackHandle* track) {
02173
02174
02175
02176
02177
02178
02179
02180
02181
02182 MSG("NtpSR",Msg::kVerbose) << "NtpSRModule::FillNtpTrackVertex" << endl;
02183
02184 NtpSRVertex& vtx = ntptrack->vtx;
02185
02186 vtx.u = track->GetVtxU();
02187 vtx.v = track->GetVtxV();
02188 vtx.x = kCos45*(vtx.u-vtx.v);
02189 vtx.y = kCos45*(vtx.u+vtx.v);
02190 vtx.z = track->GetVtxZ();
02191 vtx.t = track->GetVtxT();
02192 vtx.plane = track->GetVtxPlane();
02193 vtx.dcosu = track->GetDirCosU();
02194 vtx.dcosv = track->GetDirCosV();
02195 vtx.dcosx = kCos45*(vtx.dcosu-vtx.dcosv);
02196 vtx.dcosy = kCos45*(vtx.dcosu+vtx.dcosv);
02197 vtx.dcosz = track->GetDirCosZ();
02198
02199 const CandFitTrackHandle* fittrack
02200 = dynamic_cast<const CandFitTrackHandle*>(track);
02201 if ( fittrack ) {
02202 vtx.eu = fittrack->GetVtxUError();
02203 vtx.ev = fittrack->GetVtxVError();
02204 vtx.ex = kCos45*sqrt(vtx.eu*vtx.eu+vtx.ev*vtx.ev);
02205 vtx.ey = kCos45*sqrt(vtx.eu*vtx.eu+vtx.ev*vtx.ev);
02206 Double_t edudz = fittrack->GetVtxdUError();
02207 Double_t edvdz = fittrack->GetVtxdVError();
02208
02209
02210 vtx.edcosz = fabs(vtx.dcosz)*sqrt(pow(vtx.dcosu*edudz,2)+
02211 pow(vtx.dcosv*edvdz,2));
02212 vtx.edcosu = sqrt(fabs(vtx.dcosz))*sqrt(pow(vtx.dcosu*vtx.dcosv*edvdz,2)
02213 + pow((pow(vtx.dcosz,2)+pow(vtx.dcosv,2))*edudz,2));
02214 vtx.edcosv = sqrt(abs(vtx.dcosz))*sqrt(pow(vtx.dcosu*vtx.dcosv*edudz,2)
02215 + pow((pow(vtx.dcosz,2)+pow(vtx.dcosu,2))*edvdz,2));
02216
02217
02218 vtx.edcosx = kCos45*sqrt(vtx.edcosu*vtx.edcosu+vtx.edcosv*vtx.edcosv);
02219 vtx.edcosy = kCos45*sqrt(vtx.edcosu*vtx.edcosu+vtx.edcosv*vtx.edcosv);
02220 }
02221
02222 return;
02223
02224 }
02225
02226 void NtpSRModule::FillNtpTrackEnd(NtpSRTrack* ntptrack,
02227 const CandTrackHandle* track) {
02228
02229
02230
02231
02232
02233
02234
02235
02236
02237 MSG("NtpSR",Msg::kVerbose) << "NtpSRModule::FillNtpTrackEnd" << endl;
02238
02239 NtpSRVertex& end = ntptrack->end;
02240
02241 end.u = track->GetEndU();
02242 end.v = track->GetEndV();
02243 end.x = kCos45*(end.u-end.v);
02244 end.y = kCos45*(end.u+end.v);
02245 end.z = track->GetEndZ();
02246 end.t = track->GetEndT();
02247 end.plane = track->GetEndPlane();
02248 end.dcosu = track->GetEndDirCosU();
02249 end.dcosv = track->GetEndDirCosV();
02250 end.dcosx = kCos45*(end.dcosu-end.dcosv);
02251 end.dcosy = kCos45*(end.dcosu+end.dcosv);
02252 end.dcosz = track->GetEndDirCosZ();
02253 const CandFitTrackHandle* fittrack
02254 =dynamic_cast<const CandFitTrackHandle*>(track);
02255 if ( fittrack ) {
02256 end.eu = fittrack->GetEndUError();
02257 end.ev = fittrack->GetEndVError();
02258 end.ex = kCos45*sqrt(end.eu*end.eu+end.ev*end.ev);
02259 end.ey = kCos45*sqrt(end.eu*end.eu+end.ev*end.ev);
02260 Double_t edudz = fittrack->GetEnddUError();
02261 Double_t edvdz = fittrack->GetEnddVError();
02262
02263
02264 end.edcosz = fabs(end.dcosz)*sqrt(pow(end.dcosu*edudz,2)+
02265 pow(end.dcosv*edvdz,2));
02266
02267 end.edcosu = sqrt(fabs(end.dcosz))*sqrt(pow(end.dcosu*end.dcosv*edvdz,2)
02268 + pow((pow(end.dcosz,2)+pow(end.dcosv,2))*edudz,2));
02269 end.edcosv = sqrt(fabs(end.dcosz))*sqrt(pow(end.dcosu*end.dcosv*edudz,2)
02270 + pow((pow(end.dcosz,2)+pow(end.dcosu,2))*edvdz,2));
02271 end.edcosx = kCos45*sqrt(end.edcosu*end.edcosu+end.edcosv*end.edcosv);
02272 end.edcosy = kCos45*sqrt(end.edcosu*end.edcosu+end.edcosv*end.edcosv);
02273 }
02274
02275 return;
02276
02277 }
02278
02279 void NtpSRModule::FillNtpTrackLinearFit(NtpSRTrack* ntptrack,
02280 const CandTrackHandle* track) {
02281
02282
02283
02284
02285
02286
02287
02288
02289
02290
02291 MSG("NtpSR",Msg::kVerbose) << "NtpSRModule::FillNtpTrackLinearFit" << endl;
02292
02293
02294
02295 const Int_t nplane = 1000;
02296 Double_t uzfit[nplane]={0},vzfit[nplane]={0},ufit[nplane]={0};
02297 Double_t vfit[nplane]={0},uwfit[nplane]={0},vwfit[nplane]={0};
02298 Double_t uph[nplane]={0},vph[nplane]={0};
02299
02300
02301 TIter trackstripItr(track->GetDaughterIterator());
02302 while(CandStripHandle* trackstrip = dynamic_cast<CandStripHandle*>
02303 (trackstripItr())) {
02304
02305 Float_t tpos = trackstrip->GetTPos();
02306 Float_t phend[2] = {trackstrip->GetCharge(StripEnd::kNegative),
02307 trackstrip->GetCharge(StripEnd::kPositive)};
02308 Float_t phsum = phend[0] + phend[1];
02309
02310 Int_t iplane = trackstrip->GetPlane();
02311 assert(iplane >= 0 && iplane < nplane);
02312
02313 if ( trackstrip->GetPlaneView()==PlaneView::kU ) {
02314 uzfit[iplane] = trackstrip->GetZPos();
02315 ufit[iplane] += tpos*phsum;
02316 uph[iplane] += phsum;
02317 uwfit[iplane] = 1;
02318 }
02319 else if ( trackstrip->GetPlaneView()==PlaneView::kV ) {
02320 vzfit[iplane] = trackstrip->GetZPos();
02321 vfit[iplane] += tpos*phsum;
02322 vph[iplane] += phsum;
02323 vwfit[iplane] = 1;
02324 }
02325
02326 }
02327
02328
02329 NtpSRVertex& lin = ntptrack->lin;
02330 lin.z = track->GetVtxZ();
02331 lin.t = track->GetVtxT();
02332 Int_t timesign = 1;
02333 if ( track->GetTimeSlope() < 0. ) timesign = -1;
02334
02335 Int_t nhitplane[2] = {0};
02336 for ( Int_t ip = 0; ip < nplane; ip++ ) {
02337 if ( uph[ip] > 0 ) { nhitplane[0]++; ufit[ip] /= uph[ip]; }
02338 if ( vph[ip] > 0 ) { nhitplane[1]++; vfit[ip] /= vph[ip]; }
02339 }
02340
02341 Double_t uparm[2]={0},ueparm[2]={0},vparm[2]={0},veparm[2]={0};
02342 if ( nhitplane[0] > 1 ) {
02343 LinearFit::Weighted(nplane,uzfit,ufit,uwfit,uparm,ueparm);
02344 }
02345 if ( nhitplane[1] > 1 ) {
02346 LinearFit::Weighted(nplane,vzfit,vfit,vwfit,vparm,veparm);
02347 }
02348
02349 if ( nhitplane[0] > 1 || nhitplane[1] > 1 ) {
02350
02351 Double_t dvdz = vparm[1];
02352 lin.v = vparm[0]+lin.z*dvdz;
02353 Double_t dudz = uparm[1];
02354 lin.u = uparm[0]+lin.z*dudz;
02355 Double_t anorm = sqrt(1.+dudz*dudz+dvdz*dvdz);
02356 lin.dcosu = timesign*dudz/anorm;
02357 lin.dcosv = timesign*dvdz/anorm;
02358 lin.dcosz = timesign*1./anorm;
02359 }
02360
02361 if ( nhitplane[0] > 1 && nhitplane[1] > 1 ) {
02362
02363 Double_t dxdz = kCos45*(uparm[1]-vparm[1]);
02364 Double_t dydz = kCos45*(uparm[1]+vparm[1]);
02365 Double_t anorm = sqrt(1.+dxdz*dxdz+dydz*dydz);
02366 lin.dcosx = timesign*dxdz/anorm;
02367 lin.dcosy = timesign*dydz/anorm;
02368 lin.x = kCos45*(uparm[0]-vparm[0])+lin.z*dxdz;
02369 lin.y = kCos45*(uparm[0]+vparm[0])+lin.z*dydz;
02370 }
02371
02372 return;
02373
02374 }
02375
02376 void NtpSRModule::FillNtpSubShowerSummary(NtpSRShower *ntpshower,
02377 const CandShowerSRHandle *showerSR,
02378 const CandTrackListHandle *tracklist)
02379 {
02380 MSG("NtpSR",Msg::kVerbose) << "NtpSRModule::FillNtpSubShowerSummary" << endl;
02381
02382 Float_t totChargeU = 0, totChargeV = 0;
02383 ntpshower->sss.Zero();
02384
02385 for (int i = 0; i<showerSR->GetLastSubShower()+1; i++) {
02386 const CandSubShowerSRHandle* showercluster=showerSR->GetSubShower(i);
02387
02388 if(showercluster->GetPlaneView()==PlaneView::kU) {
02389 totChargeU += showercluster->GetEnergy();
02390 ntpshower->sss.compactU += TMath::Power(showercluster->GetEnergy(),2);
02391 }
02392 if(showercluster->GetPlaneView()==PlaneView::kV) {
02393 totChargeV += showercluster->GetEnergy();
02394 ntpshower->sss.compactV += TMath::Power(showercluster->GetEnergy(),2);
02395 }
02396
02397 if(showercluster->GetClusterID()==ClusterType::kTrkLike &&
02398 showercluster->GetNStrip()>2) {
02399 if(showercluster->GetPlaneView()==PlaneView::kU) {
02400 ntpshower->sss.nTrkLikeU += showercluster->GetNStrip();
02401 ntpshower->sss.phTrkLikeU += showercluster->GetEnergy();
02402 }
02403 else if(showercluster->GetPlaneView()==PlaneView::kV) {
02404 ntpshower->sss.nTrkLikeV += showercluster->GetNStrip();
02405 ntpshower->sss.phTrkLikeV += showercluster->GetEnergy();
02406 }
02407 }
02408 else if(showercluster->GetClusterID()==ClusterType::kEMLike) {
02409 if(showercluster->GetPlaneView()==PlaneView::kU) {
02410 ntpshower->sss.probEMU += ( showercluster->GetProbEM() *
02411 showercluster->GetEnergy() );
02412 }
02413 else if(showercluster->GetPlaneView()==PlaneView::kV) {
02414 ntpshower->sss.probEMV += ( showercluster->GetProbEM() *
02415 showercluster->GetEnergy() );
02416 }
02417 }
02418 }
02419
02420 if(totChargeU>0) {
02421 ntpshower->sss.probEMU /= totChargeU;
02422 if(ntpshower->nUcluster==1) ntpshower->sss.compactU = 1;
02423 else {
02424 ntpshower->sss.compactU /= (totChargeU*totChargeU*float(ntpshower->nUcluster));
02425 ntpshower->sss.compactU -= TMath::Power(1./float(ntpshower->nUcluster),2);
02426 if(ntpshower->sss.compactU>0)
02427 ntpshower->sss.compactU = 2.*TMath::Sqrt(ntpshower->sss.compactU);
02428 else ntpshower->sss.compactU = 0;
02429 }
02430 }
02431 if(totChargeV>0) {
02432 ntpshower->sss.probEMV /= totChargeV;
02433 if(ntpshower->nVcluster==1) ntpshower->sss.compactV = 1;
02434 else {
02435 ntpshower->sss.compactV /= (totChargeV*totChargeV*float(ntpshower->nVcluster));
02436 ntpshower->sss.compactV -= TMath::Power(1./float(ntpshower->nVcluster),2);
02437 if(ntpshower->sss.compactV>0)
02438 ntpshower->sss.compactV = 2.*TMath::Sqrt(ntpshower->sss.compactV);
02439 else ntpshower->sss.compactV = 0;
02440 }
02441 }
02442
02443 if(tracklist) {
02444 TIter trkItr(tracklist->GetDaughterIterator());
02445 for(int i=0;i<showerSR->GetLastSubShower()+1;i++) {
02446 const CandSubShowerSRHandle* showercluster = showerSR->GetSubShower(i);
02447 if(!(showercluster->GetClusterID()==ClusterType::kTrkLike &&
02448 showercluster->GetNStrip()>2)) continue;
02449 TIter clustpItr(showercluster->GetDaughterIterator());
02450 while(CandStripHandle *clustrip =
02451 dynamic_cast<CandStripHandle*>(clustpItr())){
02452 Bool_t foundStrip = false;
02453 trkItr.Reset();
02454 while(CandTrackHandle *track =
02455 dynamic_cast<CandTrackHandle*>(trkItr())){
02456 if(!track->GetCandSlice()->
02457 IsEquivalent(showerSR->GetCandSlice())) continue;
02458 TIter stpItr(track->GetDaughterIterator());
02459 while(CandStripHandle *strip =
02460 dynamic_cast<CandStripHandle*>(stpItr())){
02461 if(strip->IsEquivalent(clustrip)) {
02462 if(strip->GetPlaneView()==PlaneView::kU)
02463 ntpshower->sss.nRecoTrkU += 1;
02464 else if(strip->GetPlaneView()==PlaneView::kV)
02465 ntpshower->sss.nRecoTrkV += 1;
02466 foundStrip = true;
02467 break;
02468 }
02469 }
02470 if(foundStrip) break;
02471 }
02472 }
02473 }
02474 }
02475 }
02476
02477 void NtpSRModule::FillNtpBleach(NtpSRBleach& bleach,
02478 const CandEventHandle* cndevt,
02479 const CandRecord* cndrec,
02480 const RawRecord* rawrec) {
02481
02482
02483
02484
02485
02486
02487 bleach.lateBucketPHFraction
02488 = NtpSRBleachFiller::GetlateBucketPHFraction(cndevt,rawrec);
02489 bleach.straightPHFraction
02490 = NtpSRBleachFiller::GetstraightPHFraction(cndevt, cndrec);
02491 bleach.fixedWindowPH
02492 = NtpSRBleachFiller::GetFixedWindowPH(cndevt,cndrec);
02493 bleach.eventDuration
02494 = NtpSRBleachFiller::GetEventDuration(cndevt);
02495 return;
02496
02497 }
02498
02499 void NtpSRModule::FillNtpDataQuality(NtpSRDataQuality& ntpdataquality, TClonesArray& ntpdeadchips, const CandRecord* cndrec) {
02500
02501 MSG("NtpSR",Msg::kDebug) << "NtpSRModule::FillNtpDataQuality" << endl;
02502
02503 const CandDataQualityHandle *dqhandle
02504 = dynamic_cast <const CandDataQualityHandle*>
02505 (cndrec->FindCandHandle("CandDataQualityHandle"));
02506 if ( !dqhandle ) return;
02507
02508 ntpdataquality.trigsource = dqhandle->GetTriggerSource();
02509 ntpdataquality.trigtime = dqhandle->GetTriggerTime();
02510 ntpdataquality.errorcode = dqhandle->GetErrorCode();
02511 ntpdataquality.cratemask = dqhandle->GetCrateMask();
02512 ntpdataquality.pretrigdigits = dqhandle->GetPreTriggerDigits();
02513 ntpdataquality.posttrigdigits = dqhandle->GetPostTriggerDigits();
02514 ntpdataquality.snarlmultiplicity = dqhandle->GetSnarlMultiplicity();
02515 ntpdataquality.spillstatus = dqhandle->GetSpillStatus();
02516 ntpdataquality.spilltype = dqhandle->GetSpillType();
02517 ntpdataquality.spilltimeerror = dqhandle->GetSpillTimeError();
02518 ntpdataquality.litrigger = dqhandle->GetLiTrigger();
02519 ntpdataquality.litime = dqhandle->GetLiTime();
02520 ntpdataquality.lisubtractedtime = dqhandle->GetLiSubtractedTime();
02521 ntpdataquality.lirelativetime = dqhandle->GetLiRelativeTime();
02522 ntpdataquality.licalibpoint = dqhandle->GetLiCalibPoint();
02523 ntpdataquality.licalibtype = dqhandle->GetLiCalibType();
02524 ntpdataquality.libox = dqhandle->GetLiPulserBox();
02525 ntpdataquality.liled = dqhandle->GetLiPulserLed();
02526 ntpdataquality.lipulseheight = dqhandle->GetLiPulseHeight();
02527 ntpdataquality.lipulsewidth = dqhandle->GetLiPulseWidth();
02528 ntpdataquality.coldchips = dqhandle->GetColdChips();
02529 ntpdataquality.hotchips = dqhandle->GetHotChips();
02530 ntpdataquality.busychips = dqhandle->GetBusyChips();
02531 ntpdataquality.readouterrors = dqhandle->GetReadoutErrors();
02532 ntpdataquality.dataqualityword = (Int_t)dqhandle->GetDataQuality();
02533
02534 MSG("NtpSR",Msg::kDebug) << ntpdataquality << endl;
02535
02536 Int_t ndeadchips = 0;
02537 TIter DeadChipItr(dqhandle->GetDaughterIterator());
02538 while ( CandDeadChipHandle* deadchip=dynamic_cast<CandDeadChipHandle*>(DeadChipItr())) {
02539
02540
02541 NtpSRDeadChip* ntpdeadchip = new((ntpdeadchips)[ndeadchips++])NtpSRDeadChip();
02542 VldContext* vldc = (VldContext*)(dqhandle->GetVldContext());
02543 RawChannelId rawch = deadchip->GetChannelId();
02544 PlexHandle plex(*vldc);
02545
02546 Int_t tempchannelid=-1;
02547 Int_t plane0=-1;
02548 Int_t plane1=-1;
02549 Int_t shield=-1;
02550
02551 if( rawch.GetElecType()==ElecType::kVA
02552 && vldc->GetDetector()==Detector::kFar ){
02553
02554 Int_t crate = rawch.GetCrate();
02555 Int_t varc = rawch.GetVarcId();
02556 Int_t vmm = rawch.GetVmm();
02557 Int_t vaadc = rawch.GetVaAdcSel();
02558 Int_t vachip = rawch.GetVaChip();
02559
02560 tempchannelid=108*crate+36*varc+6*vmm+3*vaadc+vachip;
02561
02562 RawChannelId rcidCheck = RawChannelId(Detector::kFar,ElecType::kVA,
02563 crate, varc, vmm, vaadc, vachip,
02564 10);
02565
02566 if( plex.GetSEIdAltL(rcidCheck).IsValid()
02567 && plex.GetSEIdAltL(rcidCheck).GetPlane()>=0
02568 && plex.GetSEIdAltL(rcidCheck).GetPlane()<1000 ){
02569
02570 if( plex.GetSEIdAltL(rcidCheck).IsVetoShield() ){
02571 shield = plex.GetSEIdAltL(rcidCheck).GetPlane();
02572 }
02573
02574 if( !plex.GetSEIdAltL(rcidCheck).IsVetoShield() ){
02575 RawChannelId rcidLow = RawChannelId(Detector::kFar,ElecType::kVA,
02576 crate,varc,vmm,vaadc,0,2);
02577
02578 RawChannelId rcidHigh = RawChannelId(Detector::kFar,ElecType::kVA,
02579 crate,varc,vmm,vaadc,1,17);
02580
02581 if(vachip!=1) plane0 = plex.GetSEIdAltL(rcidLow).GetPlane();
02582 if(vachip!=0) plane1 = plex.GetSEIdAltL(rcidHigh).GetPlane();
02583 }
02584 }
02585
02586 }
02587
02588 if( rawch.GetElecType()==ElecType::kQIE
02589 && vldc->GetDetector()==Detector::kNear ){
02590
02591 Int_t crate = rawch.GetCrate();
02592 Int_t master = rawch.GetMaster();
02593 Int_t minder = rawch.GetMinder();
02594 Int_t menu = rawch.GetMenu();
02595
02596 plane0 = plex.GetSEIdAltL(rawch).GetPlane();
02597 tempchannelid=2560*crate+128*master+16*minder+menu;
02598 }
02599
02600
02601 ntpdeadchip->channelid = tempchannelid;
02602 ntpdeadchip->plane0 = plane0;
02603 ntpdeadchip->plane1 = plane1;
02604 ntpdeadchip->shield = shield;
02605 ntpdeadchip->errorcode = deadchip->GetErrorCode();
02606 ntpdeadchip->status = (Int_t)deadchip->GetChipStatus();
02607
02608 }
02609
02610 }
02611
02612 void NtpSRModule::FillNtpWindow(NtpSREvent* ntpevent,
02613 const CandRecord* cndrec) {
02614
02615 const VldContext& vld = *(cndrec->GetVldContext());
02616 Bool_t isND = false;
02617 if(vld.GetDetector()==Detector::kNear) isND = true;
02618
02619 ntpevent->win.begplane = ntpevent->vtx.plane - fWindowPlaneExtn;
02620 ntpevent->win.endplane = ntpevent->end.plane + fWindowPlaneExtn;
02621 ntpevent->win.begtime = ntpevent->vtx.t - fWindowTimeExtn;
02622 ntpevent->win.endtime = ntpevent->end.t + fWindowTimeExtn;
02623
02624 ntpevent->win.totalQ = 0;
02625 ntpevent->win.specQ = 0;
02626 ntpevent->win.pinstQ = 0;
02627 ntpevent->win.utotalQ = 0;
02628 ntpevent->win.uspecQ = 0;
02629 ntpevent->win.upinstQ = 0;
02630
02631 const CandStripListHandle *striplisthandle
02632 = dynamic_cast <const CandStripListHandle*>
02633 (cndrec->FindCandHandle("CandStripListHandle"));
02634 if(!striplisthandle) return;
02635
02636 TIter stripItr(striplisthandle->GetDaughterIterator());
02637 while ( CandStripHandle* strip=dynamic_cast<CandStripHandle*>(stripItr())) {
02638 Int_t uid = strip->GetUidInt();
02639 Int_t cur_plane = strip->GetPlane();
02640 Double_t cur_time = strip->GetTime();
02641
02642 if(cur_plane>=ntpevent->win.begplane &&
02643 cur_plane<=ntpevent->win.endplane &&
02644 cur_time>=ntpevent->win.begtime &&
02645 cur_time<=ntpevent->win.endtime){
02646 Int_t cur_strip = strip->GetStrip();
02647 Float_t cur_charge = strip->GetCharge(CalDigitType::kSigCorr);
02648 Bool_t isSpec = false;
02649 Bool_t isPInst = false;
02650 if(isND){
02651 if(cur_plane>FIRST_ND_SPECTROMETER_PLANE) isSpec = true;
02652 if(strip->GetPlaneView()==PlaneView::kU) {
02653 PlexPlaneId ppid(Detector::kNear,strip->GetPlane(),kFALSE);
02654 if(ppid.GetPlaneCoverage()==PlaneCoverage::kNearFull){
02655 if(cur_strip>=FIRST_FULL_PINST_STRIP_U &&
02656 cur_strip<=LAST_FULL_PINST_STRIP_U) isPInst = true;
02657 }
02658 else isPInst = true;
02659 }
02660 else if(strip->GetPlaneView()==PlaneView::kV) {
02661
02662 if(cur_strip>=FIRST_PINST_STRIP_V &&
02663 cur_strip<=LAST_PINST_STRIP_V) isPInst = true;
02664 }
02665 }
02666
02667 ntpevent->win.totalQ += cur_charge;
02668 if(isND) {
02669 if(isSpec) ntpevent->win.specQ += cur_charge;
02670 else if(isPInst) ntpevent->win.pinstQ += cur_charge;
02671 }
02672
02673 std::map<int,bool>::iterator uidItr = fUnassocStripUidMap.find(uid);
02674 if ( uidItr != fUnassocStripUidMap.end() ) {
02675 if(uidItr->second){
02676
02677 ntpevent->win.utotalQ += cur_charge;
02678 if(isND) {
02679 if(isSpec) ntpevent->win.uspecQ += cur_charge;
02680 else if(isPInst) ntpevent->win.upinstQ += cur_charge;
02681 }
02682 }
02683 }
02684 }
02685 }
02686 }
02687
02688 void NtpSRModule::FillNtpDetStatus(NtpSRDetStatus& ntpdetstatus,
02689 const VldContext& vldc) {
02690
02691
02692
02693
02694
02695
02696
02697
02698
02699
02700 MSG("NtpSR",Msg::kDebug) << "NtpSRModule::FillNtpSRDetStatus" << endl;
02701
02702
02703 ntpdetstatus.coilstatus = 0;
02704 ntpdetstatus.dcscoilstatus = CoilStatus::kUnknown;
02705 ntpdetstatus.coilcurrent1 = -999.9;
02706 ntpdetstatus.coilcurrent2 = -999.9;
02707
02708 switch (vldc.GetDetector()) {
02709
02710 case Detector::kFar:
02711 {
02712 DbiResultPtr<BfldDbiCoilState> farcoilTable(vldc);
02713 Int_t nrow = farcoilTable.GetNumRows();
02714 for ( Int_t irow = 0; irow < nrow; irow++ ) {
02715 const BfldDbiCoilState* farcoilState = farcoilTable.GetRow(irow);
02716 if ( irow==0 ) ntpdetstatus.dcscoilstatus = farcoilState->GetStatus();
02717 else {
02718
02719 if (ntpdetstatus.dcscoilstatus != (Short_t)farcoilState->GetStatus())
02720 ntpdetstatus.dcscoilstatus = CoilStatus::kUnknown;
02721 }
02722
02723 if( farcoilState->GetSupermodule()==1 )
02724 ntpdetstatus.coilcurrent1 = farcoilState->GetCurrent();
02725 else if ( farcoilState->GetSupermodule()==2 )
02726 ntpdetstatus.coilcurrent2 = farcoilState->GetCurrent();
02727 }
02728 }
02729 break;
02730
02731 case Detector::kNear:
02732 {
02733 if ( CoilTools::IsOK(vldc) )
02734 ntpdetstatus.dcscoilstatus = CoilStatus::kOK;
02735 else ntpdetstatus.dcscoilstatus = CoilStatus::kBad;
02736 if ( CoilTools::IsReverse(vldc) )
02737 ntpdetstatus.dcscoilstatus |= CoilStatus::kReverse;
02738
02739
02740 DbiResultPtr<Dcs_Mag_Near> nearcoilTable(vldc);
02741 if ( nearcoilTable.GetNumRows() > 0 ) {
02742 const Dcs_Mag_Near* nearcoilState = nearcoilTable.GetRow(0);
02743 ntpdetstatus.coilcurrent1 = nearcoilState -> GetCurrent();
02744 }
02745 }
02746 break;
02747
02748 default:
02749 break;
02750 }
02751
02752
02753
02754 ntpdetstatus.dbuhvstatus = -1;
02755 ntpdetstatus.coldchips1 = -1;
02756 ntpdetstatus.coldchips2 = -1;
02757
02758 DbiResultPtr<DbuHvFromSingles> hvTable(vldc,fHvSinglesTask);
02759 Int_t nrows = hvTable.GetNumRows();
02760
02761 for( Int_t irows = 0; irows < nrows; irows++ ){
02762 const DbuHvFromSingles* hvState = hvTable.GetRow(irows);
02763
02764
02765 if( hvState->GetSupermodule()==1
02766 && hvState->GetColdChips()>ntpdetstatus.coldchips1 ){
02767 ntpdetstatus.coldchips1 = hvState->GetColdChips();
02768 }
02769
02770
02771 if( hvState->GetSupermodule()==2
02772 && hvState->GetColdChips()>ntpdetstatus.coldchips2 ){
02773 ntpdetstatus.coldchips2 = hvState->GetColdChips();
02774 }
02775
02776
02777
02778
02779 if( ( hvState->GetStatus()==0 )
02780 || ( hvState->GetStatus()==1 && ntpdetstatus.dbuhvstatus==-1 ) ) {
02781 ntpdetstatus.dbuhvstatus=hvState->GetStatus();
02782 }
02783 }
02784
02785 MSG("NtpSR",Msg::kDebug) << ntpdetstatus << endl;
02786
02787 return;
02788
02789 }
02790
02791 void NtpSRModule::FillNtpDmxStatus(NtpSRDmxStatus& ntpdmxstatus,
02792 const CandRecord* cndrec) {
02793
02794
02795
02796
02797
02798
02799
02800
02801
02802
02803 MSG("NtpSR",Msg::kDebug) << "NtpSRModule::FillNtpDmxStatus" << endl;
02804
02805 const CandDeMuxDigitListHandle *dlh
02806 = dynamic_cast <const CandDeMuxDigitListHandle*>
02807 (cndrec -> FindCandHandle("CandDeMuxDigitListHandle"));
02808 if ( !dlh ) return;
02809
02810 Int_t demuxflagword = dlh->GetDeMuxDigitListFlagWord();
02811 if (demuxflagword & CandDeMuxDigitList::kMultipleMuonEvent)
02812 ntpdmxstatus.ismultimuon = 1;
02813
02814 if (demuxflagword & CandDeMuxDigitList::kNonPhysicalStripSolution)
02815 ntpdmxstatus.nonphysicalfail = 1;
02816
02817 if (demuxflagword & CandDeMuxDigitList::kTooFewValidPlanes)
02818 ntpdmxstatus.validplanesfail = 1;
02819
02820 if (demuxflagword & CandDeMuxDigitList::kNoVertex)
02821 ntpdmxstatus.vertexplanefail = 1;
02822
02823 ntpdmxstatus.ustrayplanes = dlh->GetNumStrayPlanesU();
02824 ntpdmxstatus.vstrayplanes = dlh->GetNumStrayPlanesV();
02825 ntpdmxstatus.uvalidplanes = dlh->GetNumValidPlanesU();
02826 ntpdmxstatus.vvalidplanes = dlh->GetNumValidPlanesV();
02827 ntpdmxstatus.avgtimeoffset = dlh->GetAvgTimeOffset();
02828
02829 MSG("NtpSR",Msg::kDebug) << ntpdmxstatus << endl;
02830
02831 return;
02832
02833 }
02834
02835 void NtpSRModule::FillNtpEventSummary(NtpSREventSummary& ntpeventsummary,
02836 const CandRecord* cndrec,
02837 const RawRecord* rawrec) {
02838
02839
02840
02841
02842
02843
02844
02845
02846
02847
02848 MSG("NtpSR",Msg::kVerbose) << "NtpSRModule::FillNtpEventSummary" << endl;
02849 const CandDigitListHandle *dlh = dynamic_cast <const CandDigitListHandle*>
02850 (cndrec->FindCandHandle("CandDigitListHandle","canddigitlist"));
02851 if ( !dlh ) return;
02852
02853 ntpeventsummary.nstrip = fStripUidMap.size();
02854 ntpeventsummary.nslice = fSliceUidMap.size();
02855 ntpeventsummary.ncluster = fClusterUidMap.size();
02856 ntpeventsummary.nshower = fShowerUidMap.size();
02857 ntpeventsummary.ntrack = fTrackUidMap.size();
02858 ntpeventsummary.nevent = fEventUidMap.size();
02859
02860 ntpeventsummary.trigtime = dlh->GetAbsTime();
02861
02862 const VldContext& vld = *(cndrec->GetVldContext());
02863 PlexHandle plex(vld);
02864
02865 Double_t maxliTime = -1;
02866 if ( rawrec ) {
02867 TIter blockIter = rawrec -> GetRawBlockIter();
02868 TObject* blockobj = 0;
02869 while ( ( blockobj = blockIter.Next() ) ) {
02870 const RawDigitDataBlock* rddb
02871 = dynamic_cast<const RawDigitDataBlock*>(blockobj);
02872 if ( rddb != 0 ) {
02873 TIter digititer = rddb -> GetDatumIter();
02874 TObject* digitobj = 0;
02875 while ( ( digitobj = digititer.Next() ) ) {
02876 RawDigit* rawdigit = dynamic_cast<RawDigit*>(digitobj);
02877 if ( rawdigit ) {
02878 RawChannelId rcid = rawdigit->GetChannel();
02879 ReadoutType::Readout_t type = plex.GetReadoutType(rcid);
02880 if ( type == ReadoutType::kFlashTrigPMT ) {
02881 if ( rawdigit->GetADC() > 100 ) {
02882 Double_t liTime =
02883 Calibrator::Instance().GetTimeFromTDC(rawdigit->GetTDC(),
02884 rawdigit->GetChannel());
02885 maxliTime = TMath::Max(liTime,maxliTime);
02886 }
02887 }
02888 }
02889 }
02890 }
02891 }
02892 }
02893 else {
02894 MSG("NtpSR",Msg::kWarning) << "Missing RawRecord!"
02895 << "\n evthdr.litime will not be properly filled." << endl;
02896 }
02897 ntpeventsummary.litime = maxliTime;
02898
02899
02900 UInt_t year,month,day,hour,minute,sec;
02901 vld.GetTimeStamp().GetDate(kTRUE,0,&year,&month,&day);
02902 vld.GetTimeStamp().GetTime(kTRUE,0,&hour,&minute,&sec);
02903 ntpeventsummary.date.year = (UShort_t)year;
02904 ntpeventsummary.date.month = (Char_t)month;
02905 ntpeventsummary.date.day = (Char_t)day;
02906 ntpeventsummary.date.hour = (Char_t)hour;
02907 ntpeventsummary.date.minute = (Char_t)minute;
02908 ntpeventsummary.date.sec = (Double_t)sec;
02909 ntpeventsummary.date.sec += ntpeventsummary.trigtime;
02910 ntpeventsummary.date.utc = vld.GetTimeStamp().GetSec();
02911
02912 Int_t minplaneall = -1;
02913 Int_t minplaneallu = -1;
02914 Int_t minplaneallv = -1;
02915 Int_t maxplaneall = -1;
02916 Int_t maxplaneallu = -1;
02917 Int_t maxplaneallv = -1;
02918 std::map<Int_t,Bool_t> planeoccupancyall;
02919 std::map<Int_t,Bool_t> planeoccupancyallu;
02920 std::map<Int_t,Bool_t> planeoccupancyallv;
02921 Float_t planepe[1000][2] = {{0},{0}};
02922
02923
02924 TIter digitItr(dlh -> GetDaughterIterator());
02925 while (CandDigitHandle* digit = dynamic_cast<CandDigitHandle*>(digitItr())) {
02926
02927
02928 PlexSEIdAltL pseid(digit->GetPlexSEIdAltL());
02929 if (pseid.IsVetoShield()) continue;
02930 ntpeventsummary.ndigit++;
02931
02932 ntpeventsummary.ph.raw += digit->GetCharge(CalDigitType::kNone);
02933 Float_t calcharge[3] = {0};
02934 if ( pseid.GetDemuxVetoFlag() == 0 ) {
02935
02936 calcharge[0] = digit->GetCharge(CalDigitType::kSigLin);
02937 calcharge[1] = digit->GetCharge(CalDigitType::kSigCorr);
02938 calcharge[2] = digit->GetCharge(CalDigitType::kPE);
02939 }
02940 else if ( pseid.GetSize() > 0 ) {
02941
02942 calcharge[0] = pseid[0].GetSigLin();
02943 calcharge[1] = pseid[0].GetSigCorr();
02944 calcharge[2] = pseid[0].GetPE();
02945 }
02946 ntpeventsummary.ph.siglin += calcharge[0];
02947 ntpeventsummary.ph.sigcor += calcharge[1];
02948 ntpeventsummary.ph.pe += calcharge[2];
02949
02950
02951
02952 Int_t iplane = pseid.GetPlane();
02953 if ( iplane < 0 || iplane >= 1000 ) continue;
02954
02955 if ( minplaneall < 0 || iplane < minplaneall ) minplaneall = iplane;
02956 if ( maxplaneall < 0 || iplane > maxplaneall ) maxplaneall = iplane;
02957 planeoccupancyall[iplane] = kTRUE;
02958
02959 switch (pseid.GetPlaneView()) {
02960 case PlaneView::kU:
02961 if (minplaneallu < 0 || iplane < minplaneallu) minplaneallu = iplane;
02962 if (maxplaneallu < 0 || iplane > maxplaneallu) maxplaneallu = iplane;
02963 planeoccupancyallu[iplane] = kTRUE;
02964 break;
02965 case PlaneView::kV:
02966 if (minplaneallv < 0 || iplane < minplaneallv) minplaneallv = iplane;
02967 if (maxplaneallv < 0 || iplane > maxplaneallv) maxplaneallv = iplane;
02968 planeoccupancyallv[iplane] = kTRUE;
02969 break;
02970 default:
02971 break;
02972 }
02973
02974
02975
02976
02977 switch (pseid.GetEnd()) {
02978 case StripEnd::kNegative:
02979 planepe[iplane][0] += calcharge[2];
02980 break;
02981 case StripEnd::kPositive:
02982 planepe[iplane][1] += calcharge[2];
02983 break;
02984 default:
02985 break;
02986 }
02987 }
02988
02989 ntpeventsummary.planeall.beg = minplaneall;
02990 ntpeventsummary.planeall.end = maxplaneall;
02991 ntpeventsummary.planeall.begu = minplaneallu;
02992 ntpeventsummary.planeall.endu = maxplaneallu;
02993 ntpeventsummary.planeall.begv = minplaneallv;
02994 ntpeventsummary.planeall.endv = maxplaneallv;
02995 std::map<Int_t,Bool_t>::iterator iter;
02996 for ( Int_t iplane = minplaneall; iplane <= maxplaneall; iplane ++ ) {
02997 iter = planeoccupancyall.find(iplane);
02998 if ( iter != planeoccupancyall.end() ) ntpeventsummary.planeall.n++;
02999 iter = planeoccupancyallu.find(iplane);
03000 if ( iter != planeoccupancyallu.end() ) ntpeventsummary.planeall.nu++;
03001 iter = planeoccupancyallv.find(iplane);
03002 if ( iter != planeoccupancyallv.end() ) ntpeventsummary.planeall.nv++;
03003 }
03004
03005
03006
03007
03008
03009 Int_t minplane = -1;
03010 Int_t maxplane = -1;
03011 bool found(0);
03012 for (Int_t iplane=minplaneall;iplane<=maxplaneall-3 && !found; iplane++) {
03013 if ( !found && (planepe[iplane][0] + planepe[iplane][1]) > 3.
03014 && (planepe[iplane+1][0] + planepe[iplane+1][1]) > 3.
03015 && (planepe[iplane+2][0] + planepe[iplane+2][1]) > 3.
03016 && (planepe[iplane+3][0] + planepe[iplane+3][1]) > 3. ) {
03017 found = 1;
03018 minplane = iplane;
03019 }
03020 }
03021 found = 0;
03022 for (Int_t iplane=maxplaneall;iplane>=minplaneall+3 && !found; iplane--) {
03023 if ( !found && (planepe[iplane][0] + planepe[iplane][1]) > 3.
03024 && (planepe[iplane-1][0] + planepe[iplane-1][1]) > 3.
03025 && (planepe[iplane-2][0] + planepe[iplane-2][1]) > 3.
03026 && (planepe[iplane-3][0] + planepe[iplane-3][1]) > 3. ) {
03027 found = 1;
03028 maxplane = iplane;
03029 }
03030 }
03031
03032 Int_t minplaneu = -1;
03033 Int_t minplanev = -1;
03034 Int_t maxplaneu = -1;
03035 Int_t maxplanev = -1;
03036 for (Int_t iplane = minplane; iplane <= maxplane; iplane++ ) {
03037 if ( iplane >= 0 ) {
03038 if ( planepe[iplane][0] + planepe[iplane][1] > 0. ) {
03039 ntpeventsummary.plane.n++;
03040 iter = planeoccupancyallu.find(iplane);
03041 if ( iter != planeoccupancyallu.end() ) {
03042 ntpeventsummary.plane.nu++;
03043 if ( minplaneu < 0 || iplane < minplaneu ) minplaneu = iplane;
03044 if ( maxplaneu < 0 || iplane > maxplaneu ) maxplaneu = iplane;
03045 }
03046 iter = planeoccupancyallv.find(iplane);
03047 if ( iter != planeoccupancyallv.end() ) {
03048 ntpeventsummary.plane.nv++;
03049 if ( minplanev < 0 || iplane < minplanev ) minplanev = iplane;
03050 if ( maxplanev < 0 || iplane > maxplanev ) maxplanev = iplane;
03051 }
03052 }
03053 }
03054 }
03055
03056 ntpeventsummary.plane.beg = minplane;
03057 ntpeventsummary.plane.end = maxplane;
03058 ntpeventsummary.plane.begu = minplaneu;
03059 ntpeventsummary.plane.endu = maxplaneu;
03060 ntpeventsummary.plane.begv = minplanev;
03061 ntpeventsummary.plane.endv = maxplanev;
03062
03063 MSG("NtpSR",Msg::kDebug) << ntpeventsummary << endl;
03064
03065 return;
03066
03067 }
03068
03069
03070
03071
03072
03073
