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