00001
00002
00003
00004
00005
00006
00008
00009 #include "AnalysisNtuples/Module/ANtpRecoNtpManipulator.h"
00010 #include "AnalysisNtuples/Module/ANtpInfoObjectFillerNC.h"
00011 #include "AnalysisNtuples/ANtpTrackInfoNC.h"
00012 #include "AnalysisNtuples/ANtpShowerInfoNC.h"
00013 #include "AnalysisNtuples/ANtpEventInfoNC.h"
00014 #include "MessageService/MsgService.h"
00015 #include "AnalysisNtuples/ANtpTruthInfoBeam.h"
00016 #include "AnalysisNtuples/ANtpBeamInfo.h"
00017 #include "AnalysisNtuples/ANtpDefaultValue.h"
00018 #include "StandardNtuple/NtpStRecord.h"
00019 #include "CandNtupleSR/NtpSRRecord.h"
00020 #include "CandNtupleSR/NtpSREvent.h"
00021 #include "CandNtupleSR/NtpSRStrip.h"
00022 #include "CandNtupleSR/NtpSRCluster.h"
00023 #include "CandNtupleSR/NtpSRTrack.h"
00024 #include "CandNtupleSR/NtpSRShower.h"
00025 #include "MCNtuple/NtpMCTruth.h"
00026 #include "MCNtuple/NtpMCStdHep.h"
00027 #include "TruthHelperNtuple/NtpTHEvent.h"
00028 #include "DatabaseInterface/DbiResultPtr.h"
00029 #include "BField/BfieldCoilCurrent.h"
00030 #include "BeamDataUtil/BDSpillAccessor.h"
00031 #include "BeamDataUtil/BeamMonSpill.h"
00032 #include "BeamDataUtil/BMSpillAna.h"
00033 #include "BeamDataNtuple/NtpBDLiteRecord.h"
00034 #include "SpillTiming/SpillTimeFinder.h"
00035 #include "VertexFinder/NtpVtxFinder/NtpVtxFinder.h"
00036
00037
00038
00039
00040 #include "Conventions/Munits.h"
00041 #include "Conventions/PlaneView.h"
00042 #include "Conventions/ReleaseType.h"
00043
00044 #include "Registry/Registry.h"
00045
00046
00047 #include "TFile.h"
00048 #include "TH1.h"
00049 #include "TSystem.h"
00050
00051
00052 #include <cassert>
00053 #include <algorithm>
00054 #include <numeric>
00055 #include <cmath>
00056 #include <map>
00057 #include <set>
00058
00059 static int kNearPartial = 1;
00060 static int kNearFull = 2;
00061 static int kComplete = 3;
00062 static int kU = 2;
00063 static int kV = 3;
00064
00065 ClassImp(ANtpInfoObjectFillerNC)
00066
00067 CVSID("$Id: ANtpInfoObjectFillerNC.cxx,v 1.50 2010/07/28 20:51:36 rhatcher Exp $");
00068
00069
00070 ANtpInfoObjectFillerNC::ANtpInfoObjectFillerNC() :
00071
00072 ANtpInfoObjectFillerBeam(Detector::kNear),
00073 fClusterArray(0),
00074 fkNNSet(false),
00075 fVHSPlanes(20),
00076 fVHSStrips(20)
00077 {
00078 MSG("ANtpInfoObjectFillerNC", Msg::kDebug)
00079 << "ANtpInfoObjectFillerNC::Default Constructor" << endl;
00080
00081 FillStripToPixelMaps();
00082
00083 }
00084
00085
00086 ANtpInfoObjectFillerNC::ANtpInfoObjectFillerNC(Detector::Detector_t detector) :
00087 ANtpInfoObjectFillerBeam(detector),
00088 fClusterArray(0),
00089 fkNNSet(false),
00090 fVHSPlanes(20),
00091 fVHSStrips(20)
00092 {
00093 MSG("ANtpInfoObjectFillerNC", Msg::kDebug)
00094 << "ANtpInfoObjectFillerNC::Constructor" << endl;
00095
00096 FillStripToPixelMaps();
00097
00098 }
00099
00100
00101 ANtpInfoObjectFillerNC::~ANtpInfoObjectFillerNC()
00102 {
00103 if(fAnpInterfaceRO) delete fAnpInterfaceRO;
00104 if(fAnpInterfaceJM) delete fAnpInterfaceJM;
00105 }
00106
00107
00108 void ANtpInfoObjectFillerNC::FillStripToPixelMaps()
00109 {
00110
00111
00112 if(fDetector == Detector::kFar){
00113
00114
00115 Int_t stripToPixelWest[] = {0, 2, 5, 7, 8, 10, 13, 15, 0, 2, 5, 7, 1, 3, 4, 6, 9, 11, 12, 14, 1, 3, 4, 6,
00116 0, 2, 5, 7, 8, 10, 13, 15, 0, 2, 5, 7, 1, 3, 4, 6, 9, 11, 12, 14, 1, 3, 4, 6,
00117 0, 2, 5, 7, 8, 10, 13, 15, 0, 2, 5, 7, 1, 3, 4, 6, 9, 11, 12, 14, 1, 3, 4, 6,
00118 0, 2, 5, 7, 8, 10, 13, 15, 0, 2, 5, 7, 1, 3, 4, 6, 9, 11, 12, 14, 1, 3, 4, 6,
00119 0, 2, 5, 7, 8, 10, 13, 15, 0, 2, 5, 7, 1, 3, 4, 6, 9, 11, 12, 14, 1, 3, 4, 6,
00120 0, 2, 5, 7, 8, 10, 13, 15, 0, 2, 5, 7, 1, 3, 4, 6, 9, 11, 12, 14, 1, 3, 4, 6,
00121 0, 2, 5, 7, 8, 10, 13, 15, 0, 2, 5, 7, 1, 3, 4, 6, 9, 11, 12, 14, 1, 3, 4, 6,
00122 0, 2, 5, 7, 8, 10, 13, 15, 0, 2, 5, 7, 1, 3, 4, 6, 9, 11, 12, 14, 1, 3, 4, 6};
00123
00124 Int_t stripToPixelEast[] = {0, 2, 5, 7, 8, 10, 13, 15, 0, 2, 5, 7, 1, 3, 4, 6, 9, 11, 12, 14, 1, 3, 4, 6,
00125 2, 5, 7, 8, 10, 13, 15, 0, 2, 5, 7, 1, 3, 4, 6, 9, 11, 12, 14, 1, 3, 4, 6, 0,
00126 5, 7, 8, 10, 13, 15, 0, 2, 5, 7, 1, 3, 4, 6, 9, 11, 12, 14, 1, 3, 4, 6, 0, 2,
00127 7, 8, 10, 13, 15, 0, 2, 5, 7, 1, 3, 4, 6, 9, 11, 12, 14, 1, 3, 4, 6, 0, 2, 5,
00128 8, 10, 13, 15, 0, 2, 5, 7, 1, 3, 4, 6, 9, 11, 12, 14, 1, 3, 4, 6, 0, 2, 5, 7,
00129 10, 13, 15, 0, 2, 5, 7, 1, 3, 4, 6, 9, 11, 12, 14, 1, 3, 4, 6, 0, 2, 5, 7, 8,
00130 13, 15, 0, 2, 5, 7, 1, 3, 4, 6, 9, 11, 12, 14, 1, 3, 4, 6, 0, 2, 5, 7, 8, 10,
00131 15, 0, 2, 5, 7, 1, 3, 4, 6, 9, 11, 12, 14, 1, 3, 4, 6, 0, 2, 5, 7, 8, 10, 13};
00132
00133
00134 for(Int_t i = 0; i < 192; ++i){
00135 fStripToPixelMapWest[i] = stripToPixelWest[i];
00136 fStripToPixelMapEast[i] = stripToPixelEast[i];
00137 }
00138
00139
00140 for(Int_t i = 0; i < 486; ++i) fPlaneCoverage[i] = kComplete;
00141 }
00142
00143 else if(fDetector == Detector::kNear){
00144
00145
00146 Int_t pixelToStripU1[] = {40, 53, 22, 35, 11, 41, 54, 23,
00147 36, 10, 42, 55, 24, 37, 9, 43,
00148 56, 25, 38, 8, 44, 57, 26, 39,
00149 7, 45, 58, 27, 19, 6, 46, 59,
00150 28, 18, 5, 47, 60, 29, 17, 4,
00151 48, 61, 30, 16, 3, 49, 62, 31,
00152 15, 2, 50, 63, 32, 14, 1, 51,
00153 20, 33, 13, 0, 52, 21, 34, 12};
00154
00155 Int_t pixelToStripV1[] = {27, 14, 45, 32, 56, 26, 13, 44,
00156 31, 57, 25, 12, 43, 30, 58, 24,
00157 11, 42, 29, 59, 23, 10, 41, 28,
00158 60, 22, 9, 40, 48, 61, 21, 8,
00159 39, 49, 62, 20, 7, 38, 50, 63,
00160 19, 6, 37, 51, 64, 18, 5, 36,
00161 52, 65, 17, 4, 35, 53, 66, 16,
00162 47, 34, 54, 67, 15, 46, 33, 55};
00163
00164 Int_t pixelToStripU2[] = {68, 81, 94, 59, -1, 69, 82, 95,
00165 60, -1, 70, 83, 48, 61, -1, 71,
00166 84, 49, 62, -1, 72, 85, 50, 63,
00167 -1, 73, 86, 51, 64, -1, 74, 87,
00168 52, 65, -1, 75, 88, 53, 66, -1,
00169 76, 89, 54, 67, -1, 77, 90, 55,
00170 -1, -1, 78, 91, 56, -1, -1, 79,
00171 92, 57, -1, -1, 80, 93, 58, -1};
00172
00173 Int_t pixelToStripV2[] = {27, 14, 1, 36, -1, 26, 13, 0,
00174 35, -1, 25, 12, 47, 34, -1, 24,
00175 11, 46, 33, -1, 23, 10, 45, 32,
00176 -1, 22, 9, 44, 31, -1, 21, 8,
00177 43, 30, -1, 20, 7, 42, 29, -1,
00178 19, 6, 41, 28, -1, 18, 5, 40,
00179 -1, -1, 17, 4, 39, -1, -1, 16,
00180 3, 38, -1, -1, 15, 2, 37, -1};
00181
00182 Int_t pixelToStripU3[] = {47, 34, 21, 8, -1, 46, 33, 20,
00183 7, -1, 45, 32, 19, 6, -1, 44,
00184 31, 18, 5, -1, 43, 30, 17, 4,
00185 -1, 42, 29, 16, 3, -1, 41, 28,
00186 15, 2, -1, 40, 27, 14, 1, -1,
00187 39, 26, 13, 0, -1, 38, 25, 12,
00188 -1, -1, 37, 24, 11, -1, -1, 36,
00189 23, 10, -1, -1, 35, 22, 9, -1};
00190
00191 Int_t pixelToStripV3[] = {48, 61, 74, 87, -1, 49, 62, 75,
00192 88, -1, 50, 63, 76, 89, -1, 51,
00193 64, 77, 90, -1, 52, 65, 78, 91,
00194 -1, 53, 66, 79, 92, -1, 54, 67,
00195 80, 93, -1, 55, 68, 81, 94, -1,
00196 56, 69, 82, 95, -1, 57, 70, 83,
00197 -1, -1, 58, 71, 84, -1, -1, 59,
00198 72, 85, -1, -1, 60, 73, 86, -1};
00199
00200
00201 for(Int_t i = 0; i < 64; ++i){
00202 fStripToPixelMapNearU1[pixelToStripU1[i]] = i;
00203 fStripToPixelMapNearV1[pixelToStripV1[i]] = i;
00204
00205 fStripToPixelMapNearU2[i] = -1;
00206 if(pixelToStripU2[i]>-1) fStripToPixelMapNearU2[pixelToStripU2[i]] = i;
00207
00208 fStripToPixelMapNearV2[i] = -1;
00209 if(pixelToStripV2[i]>-1) fStripToPixelMapNearV2[pixelToStripV2[i]] = i;
00210
00211 fStripToPixelMapNearU3[i] = -1;
00212 if(pixelToStripU3[i]>-1) fStripToPixelMapNearU3[pixelToStripU3[i]] = i;
00213
00214 fStripToPixelMapNearV3[i] = -1;
00215 if(pixelToStripV3[i]>-1) fStripToPixelMapNearV3[pixelToStripV3[i]] = i;
00216
00217 }
00218
00219
00220 for(Int_t i = 0; i < 282; ++i){
00221 if( i%5 == 1 ) fPlaneCoverage[i] = kNearFull;
00222 else fPlaneCoverage[i] = kNearPartial;
00223 }
00224 }
00225
00226 MSG("ANtpInfoObjectFillerNC", Msg::kDebug)
00227 << "ANtpInfoObjectFillerNC::Constructor" << endl;
00228
00229 }
00230
00231
00232 void ANtpInfoObjectFillerNC::SetDetector(Detector::Detector_t detector)
00233 {
00234
00235 fDetector = detector;
00236 FillStripToPixelMaps();
00237 return;
00238 }
00239
00240
00241 void ANtpInfoObjectFillerNC::SetClusterArray(const TClonesArray *clusters) {
00242 fClusterArray = clusters;
00243 }
00244
00245
00246 bool ANtpInfoObjectFillerNC::FillInformation(int event,
00247 ANtpRecoNtpManipulator *ntpManipulator,
00248 ANtpEventInfoNC *eventInfo,
00249 ANtpTrackInfoNC *trackInfo,
00250 ANtpShowerInfoNC *showerInfo,
00251 ANtpTruthInfoBeam *truthInfo)
00252 {
00253
00254 NtpSRTrack *ntpTrack = 0;
00255 NtpSREvent *ntpEvent = 0;
00256 NtpSRShower *ntpShower = 0;
00257 NtpMCTruth *ntpMCTruth = 0;
00258 NtpTHEvent *ntpTHEvent = 0;
00259
00260 MAXMSG("ANtpInfoObjectFillerNC", Msg::kDebug,20)
00261 << "Begin Fill Information method..." << endl;
00262
00263
00264
00265
00266 ntpManipulator->GetEventManipulator()->SetEventInSnarl(event);
00267 ntpEvent = ntpManipulator->GetEventManipulator()->GetEvent();
00268 if(!ntpEvent) return false;
00269
00270 FillEventInformation(ntpManipulator, ntpEvent, eventInfo);
00271 FillEventTimingAndActivityInformation(ntpManipulator, event, eventInfo);
00272
00273
00274
00275 ntpShower = ntpManipulator->GetEventManipulator()->GetPrimaryShower();
00276 if(ntpShower) FillShowerInformation(ntpShower, showerInfo, eventInfo);
00277
00278
00279
00280 ntpTrack = ntpManipulator->GetEventManipulator()->GetPrimaryTrack();
00281 if(ntpTrack)
00282 FillTrackInformation(ntpManipulator, ntpTrack, ntpEvent, trackInfo, eventInfo);
00283
00284 FillCrossOverInformation(ntpTrack, ntpShower, ntpEvent, ntpManipulator,
00285 trackInfo, showerInfo, eventInfo);
00286
00287
00288
00289 ntpTHEvent = ntpManipulator->GetMCManipulator()->GetNtpTHEvent(event);
00290
00291 if(ntpTHEvent)
00292 ntpMCTruth = ntpManipulator->GetMCManipulator()->GetNtpMCTruth(ntpTHEvent->neumc);
00293 if(ntpMCTruth)
00294 FillMCInformation(ntpMCTruth, ntpManipulator->GetStdHepArray(), truthInfo);
00295
00296 if (truthInfo && ntpTHEvent)
00297 truthInfo->eventCompleteness = ntpTHEvent->completeall;
00298
00299 return true;
00300 }
00301
00302
00303 void ANtpInfoObjectFillerNC::
00304 GetEvtVtxWithFixup(ANtpRecoNtpManipulator* ntpManipulator,
00305 int& vtxPlane, float& vtxZ) const
00306 {
00307 ANtpEventManipulator* evtManip = ntpManipulator->GetEventManipulator();
00308 const NtpSREvent* ntpEvent = evtManip->GetEvent();
00309
00310 vtxPlane = ntpEvent->vtx.plane;
00311 vtxZ = ntpEvent->vtx.z;
00312
00313
00314 const TString reco = ntpManipulator->GetReleaseMCType();
00315 const bool isCedar = reco.Contains("Cedar");
00316 if(isCedar){
00317 NtpVtxFinder vtxf;
00318 vtxf.SetTargetEvent(ntpEvent->index, ntpManipulator->GetNtpStRecord());
00319 if(vtxf.FindVertex() > 0){
00320 vtxPlane = vtxf.VtxPlane();
00321 vtxZ = vtxf.VtxZ();
00322 }
00323 }
00324
00325 const NtpSRShower* ntpShower = evtManip->GetPrimaryShower();
00326 const NtpSRTrack* ntpTrack = evtManip->GetPrimaryTrack();
00327
00328 if(ntpTrack && ntpShower){
00329 if(ntpEvent->ntrack > 0 &&
00330 ntpTrack->plane.n > ntpShower->plane.n){
00331 vtxPlane = ntpTrack->vtx.plane;
00332 vtxZ = ntpTrack->vtx.z - 3.92 * Munits::cm;
00333 }
00334 }
00335 }
00336
00337
00338 Int_t ANtpInfoObjectFillerNC::contiguousPlanes(const NtpSREvent* ntpEvent,
00339 Float_t phPlaneCut)
00340 {
00341 TH1F* lenepl=new TH1F("lenepl","lenepl",120,0.,120.);
00342 for(int it = 0; it < ntpEvent->nstrip; ++it){
00343 const NtpSRStrip* ntpStrip = (NtpSRStrip*)fStripArray->At(ntpEvent->stp[it]);
00344 Float_t stripPH=ntpStrip->ph0.pe+ntpStrip->ph1.pe;
00345 Float_t deltaplanes = ntpStrip->plane-ntpEvent->vtx.plane;
00346
00347 if(deltaplanes>=0.)
00348 lenepl->Fill(deltaplanes-0.5,stripPH);
00349 }
00350
00351 bool PlaneCut=true;
00352 Int_t _contPlaneCount=0;
00353 if(lenepl->GetBinContent(1)>phPlaneCut) _contPlaneCount++;
00354
00355 for(int i=2;i<lenepl->GetNbinsX(); i++){
00356 if(lenepl->GetBinContent(i)>phPlaneCut && PlaneCut) _contPlaneCount++;
00357 else PlaneCut=false;
00358 }
00359 delete lenepl;
00360 return _contPlaneCount;
00361 }
00362
00363
00364 void ANtpInfoObjectFillerNC::
00365 GetStripEventTime(const NtpSREvent* ntpEvent,
00366 int vtxPlane, float vtxZ,
00367 double& evtTime,
00368 double& ToFCorrectedEvtTime) const
00369 {
00370 vector<double> stpTimes;
00371
00372 ToFCorrectedEvtTime = 0;
00373 Int_t tfcount=0;
00374
00375 for(int ns = 0; ns < ntpEvent->nstrip; ++ns){
00376 const int stpIdx = ntpEvent->stp[ns];
00377 if(stpIdx < 0) continue;
00378 const NtpSRStrip* ntpStrip = (NtpSRStrip*)fStripArray->At(stpIdx);
00379
00380 if(ntpStrip->ph1.pe >=2) {
00381 ToFCorrectedEvtTime += ntpStrip->time1-
00382 (fabs(ntpStrip->z - vtxZ)/Munits::c_light);
00383 tfcount++;
00384 }
00385
00386 if(ntpStrip->plane >= vtxPlane &&
00387 ntpStrip->plane < vtxPlane + 5){
00388 stpTimes.push_back(ntpStrip->time1);
00389 }
00390 }
00391
00392 if(tfcount) ToFCorrectedEvtTime /= (double)tfcount;
00393 else ToFCorrectedEvtTime=ANtpDefVal::kDouble;
00394
00395
00396 sort(stpTimes.begin(), stpTimes.end());
00397
00398 if(stpTimes.size() % 2)
00399 evtTime = *(stpTimes.begin()+stpTimes.size()/2);
00400 else
00401 evtTime = ( *(stpTimes.begin()+stpTimes.size()/2)
00402 + *(stpTimes.begin()+stpTimes.size()/2-1) )/2.;
00403 }
00404
00405
00406
00407 void ANtpInfoObjectFillerNC::
00408 FillEventTimingAndActivityInformation(ANtpRecoNtpManipulator* ntpManipulator,
00409 int event,
00410 ANtpEventInfoNC* eventInfo)
00411 {
00412 ANtpEventManipulator* evtManip = ntpManipulator->GetEventManipulator();
00413
00414
00415 evtManip->SetEventInSnarl(event);
00416
00417 const NtpSREvent* ntpEvent = evtManip->GetEvent();
00418
00419
00420
00421 int thisVtxPlane; float thisVtxZ;
00422 GetEvtVtxWithFixup(ntpManipulator, thisVtxPlane, thisVtxZ);
00423
00424 double thisEvtTime; double thisToFEventTime;
00425 GetStripEventTime(ntpEvent, thisVtxPlane, thisVtxZ,
00426 thisEvtTime, thisToFEventTime);
00427
00428
00429 Double_t timeDiff = ANtpDefVal::kFloat;
00430
00431 Int_t closeEventIndex = ANtpDefVal::kInt;
00432
00433 Float_t minDeltaZ = ANtpDefVal::kFloat;
00434
00435 Double_t mintime = ANtpDefVal::kFloat;
00436
00437
00438 for(int i = 0; i <= ntpManipulator->GetEventArray()->GetLast(); ++i){
00439
00440 if(i == event) continue;
00441
00442 evtManip->SetEventInSnarl(i);
00443 const NtpSREvent* otherEvent = evtManip->GetEvent();
00444
00445 int otherVtxPlane; float otherVtxZ;
00446 GetEvtVtxWithFixup(ntpManipulator, otherVtxPlane, otherVtxZ);
00447
00448 double otherEvtTime; double otherToFEventTime;
00449 GetStripEventTime(otherEvent, otherVtxPlane, otherVtxZ,
00450 otherEvtTime, otherToFEventTime);
00451
00452 if( otherToFEventTime !=ANtpDefVal::kDouble)
00453 otherToFEventTime -= ((otherVtxZ-thisVtxZ)/Munits::c_light);
00454
00455 const double deltaT = otherEvtTime - thisEvtTime;
00456 if(TMath::Abs(deltaT) < TMath::Abs(timeDiff)){
00457 timeDiff = deltaT;
00458 closeEventIndex = i;
00459 eventInfo->closeTimeDeltaZ = thisVtxZ - otherVtxZ;
00460 }
00461
00462
00463 const float deltaZ = thisVtxZ - otherVtxZ;
00464 if(TMath::Abs(deltaZ) < TMath::Abs(minDeltaZ)){
00465 minDeltaZ = deltaZ;
00466 }
00467
00468 if(thisToFEventTime != ANtpDefVal::kDouble &&
00469 otherToFEventTime != ANtpDefVal::kDouble)
00470 if(fabs(thisToFEventTime-otherToFEventTime) < fabs(mintime))
00471 mintime = thisToFEventTime-otherToFEventTime;
00472 }
00473
00474 eventInfo->minTimeSeparation = timeDiff;
00475 eventInfo->closeTimeEvent = closeEventIndex;
00476 eventInfo->minDeltaZ = minDeltaZ;
00477 eventInfo->medianTime = thisEvtTime;
00478 eventInfo->evtTimeDiff = mintime;
00479
00480
00481
00482
00483 map<int, set<int> > mplanes;
00484
00485 for(int i = 0; i <= ntpManipulator->GetEventArray()->GetLast(); ++i){
00486
00487 if(i == event) continue;
00488 evtManip->SetEventInSnarl(i);
00489 const NtpSREvent* otherEvent = evtManip->GetEvent();
00490
00491
00492 for(Int_t ns = 0; ns < otherEvent->nstrip; ++ns){
00493 const int idx = otherEvent->stp[ns];
00494 if(idx < 0) continue;
00495 NtpSRStrip* ntpStrip = (NtpSRStrip*)fStripArray->At(idx);
00496
00497 if(ntpStrip->ph1.pe >= 2)
00498 mplanes[ntpStrip->plane].insert(ntpStrip->strip);
00499 }
00500 }
00501
00502
00503 int nstrips =0;
00504
00505 int sharedst=0;
00506
00507
00508 for(int ns = 0; ns < ntpEvent->nstrip; ++ns){
00509 const int stripIdx = ntpEvent->stp[ns];
00510 if(stripIdx < 0) continue;
00511 NtpSRStrip* ntpStrip = (NtpSRStrip*)fStripArray->At(stripIdx);
00512
00513 if(ntpStrip->ph1.pe < 2.) continue;
00514 ++nstrips;
00515 map<int, set<int> >::const_iterator PlanesIter =
00516 mplanes.find(ntpStrip->plane);
00517
00518
00519 if(PlanesIter == mplanes.end()) continue;
00520 set<int> stps = PlanesIter->second;
00521 const unsigned int stripn = ntpStrip->strip;
00522
00523
00524 for(set<int>::iterator stpit = stps.begin();
00525 stpit != stps.end();
00526 ++stpit){
00527
00528 if(abs(int(stripn)-int(*stpit)) <= 1){
00529 ++sharedst;
00530 break;
00531 }
00532 }
00533 }
00534
00535 if(nstrips)
00536 eventInfo->sharedStripFraction = float(sharedst)/nstrips;
00537
00538
00539
00540 ntpManipulator->GetSliceManipulator()->SetSliceInSnarl(ntpEvent->slc);
00541 NtpSRSlice* ntpSlice = ntpManipulator->GetSliceManipulator()->GetSlice();
00542
00543 float slicePH= 0;
00544 float eventPH = 0;
00545 for(int it = 0; it < ntpSlice->nstrip; ++it){
00546 const NtpSRStrip* ntpStrip = (NtpSRStrip*)fStripArray->At(ntpSlice->stp[it]);
00547 if(ntpStrip->ph1.pe < 2) continue;
00548 slicePH += ntpStrip->ph1.sigcor;
00549 }
00550
00551 for(int it = 0; it < ntpEvent->nstrip; ++it){
00552 const NtpSRStrip* ntpStrip = (NtpSRStrip*)fStripArray->At(ntpEvent->stp[it]);
00553 if(ntpStrip->ph1.pe < 2) continue;
00554 eventPH += ntpStrip->ph1.sigcor;
00555 }
00556
00557 if(slicePH)
00558 eventInfo->slicePHFraction = eventPH/slicePH;
00559
00560
00561 Int_t consecutive=1;
00562 Int_t maxconsecutive=1;
00563
00564 vector <int> stripplanes;
00565 for(int it = 0; it < ntpEvent->nstrip; ++it){
00566 const NtpSRStrip* ntpStrip = (NtpSRStrip*)fStripArray->At(ntpEvent->stp[it]);
00567 if(ntpStrip->ph1.pe<2.) continue;
00568
00569 bool flag=false;
00570 for(unsigned int iv=0;iv<stripplanes.size();iv++)
00571 if(stripplanes[iv]==ntpStrip->plane) flag=true;
00572 if(!flag) stripplanes.push_back(ntpStrip->plane);
00573 }
00574
00575 for(unsigned int iv=0;iv<stripplanes.size();iv++)
00576 sort(stripplanes.begin(),stripplanes.end());
00577
00578 for(unsigned int iv=1;iv<stripplanes.size();iv++)
00579 if((-stripplanes[0]+stripplanes[iv]-(consecutive-1))==1)
00580 consecutive++;
00581 vector <int> max_c;
00582
00583 for(unsigned int iv=1;iv<stripplanes.size();iv++){
00584 if((stripplanes[iv]-stripplanes[iv-1])==1){
00585 maxconsecutive++;
00586 }else{
00587 max_c.push_back(maxconsecutive);
00588 maxconsecutive=1;
00589 }
00590 }
00591
00592 max_c.push_back(maxconsecutive);
00593
00594 int max=0;
00595 for(unsigned int i=0;i<max_c.size();i++){
00596 if(max_c[i]>max){
00597 max=max_c[i];
00598 }
00599 }
00600
00601 eventInfo->consecutivePlanes=consecutive;
00602 eventInfo->maxConsecutivePlanes=max;
00603
00604
00605
00606
00607
00608
00609
00610
00611
00612 eventInfo->contPlaneCount=contiguousPlanes(ntpEvent,0);
00613
00614 eventInfo->contPlaneCount050=contiguousPlanes(ntpEvent,2.0);
00615
00616
00617
00618 const double activityTimeStart = thisEvtTime - 4e-8;
00619 const double activityTimeStop = thisEvtTime + 4e-8;
00620
00621 Int_t edgeActivityStrips=0;
00622 Float_t edgeActivityPH=0;
00623 Int_t oppEdgeStrips=0;
00624 Float_t oppEdgePH=0;
00625
00626
00627 for(int i = 0; i <= fStripArray->GetLast(); ++i){
00628 const NtpSRStrip* strip = (NtpSRStrip*)fStripArray->At(i);
00629
00630
00631 if(strip->time1 <= activityTimeStart || strip->time1 >= activityTimeStop)
00632 continue;
00633
00634
00635 if(strip->plane >= 121) continue;
00636
00637 if(strip->ph1.pe < 2) continue;
00638
00639
00640 if(strip->plane % 2 == 1 && strip->tpos < -0.24){
00641 ++edgeActivityStrips;
00642 edgeActivityPH += strip->ph1.sigcor;
00643 }
00644
00645 if(strip->plane % 2 == 1 && strip->tpos > 2.27){
00646 ++oppEdgeStrips;
00647 oppEdgePH += strip->ph1.sigcor;
00648 }
00649
00650 if(strip->plane % 2 == 0 && strip->tpos > 0.24){
00651 ++edgeActivityStrips;
00652 edgeActivityPH += strip->ph1.sigcor;
00653 }
00654
00655 if(strip->plane % 2 == 0 && strip->tpos <- 2.27){
00656 ++oppEdgeStrips;
00657 oppEdgePH += strip->ph1.sigcor;
00658 }
00659 }
00660
00661 eventInfo->edgeActivityPH = edgeActivityPH;
00662 eventInfo->edgeActivityStrips = edgeActivityStrips;
00663 eventInfo->oppEdgePH = oppEdgePH;
00664 eventInfo->oppEdgeStrips = oppEdgeStrips;
00665
00666
00667
00668 int largestIndex = 0;
00669 double largestPH = 0;
00670 for(int i = 0; i <= ntpManipulator->GetEventArray()->GetLast(); ++i){
00671 evtManip->SetEventInSnarl(i);
00672 const NtpSREvent* otherEvent = evtManip->GetEvent();
00673 if(otherEvent->ph.sigcor > largestPH){
00674 largestPH = otherEvent->ph.sigcor;
00675 largestIndex = i;
00676 }
00677 }
00678
00679 if(eventInfo->index == largestIndex)
00680 eventInfo->largestEventInSnarl = 1;
00681
00682
00683 ntpManipulator->GetEventManipulator()->SetEventInSnarl(event);
00684 }
00685
00686
00687
00688 void ANtpInfoObjectFillerNC::InitializekNN(ANtpRecoNtpManipulator *ntpManipulator)
00689 {
00690
00691 if(!fkNNSet){
00692
00693
00694 TString baseConf = getenv("SRT_PUBLIC_CONTEXT");
00695
00696
00697 if(baseConf == ""){
00698 MSG("ANtpInfoObjectFillerNC", Msg::kFatal) << "no SRT_PUBLIC_CONTEXT set"
00699 << endl;
00700 assert(false);
00701 }
00702
00703 Registry *regRO = new Registry(false);
00704 regRO->Set("InterfaceConfigPath", baseConf+"/PhysicsNtuple/Config/Config2008Test.txt");
00705 regRO->Set("FillkNNFilePath", baseConf+"/NuMuBar/data/knn.physics.far.daikon_04.dogwood1.L010z185i.root");
00706 if(fDetector == Detector::kNear)
00707 regRO->Set("FillkNNFilePath", baseConf+"/NuMuBar/data/knn.physics.near.daikon_04.dogwood1.L010z185i.root");
00708
00709
00710 Registry *regJM = new Registry(false);
00711 Registry *kreg = new Registry(false);
00712
00713
00714 kreg->Set("AlgSnarlName", "FillkNN");
00715 kreg->Set("FillkNNKeySignal", int(1));
00716
00717 kreg->Set("FillkNNKNeighbor", 99);
00718 kreg->Set("FillkNNKNeighborMod", int(5));
00719 kreg->Set("FillkNNAddEvent", "no");
00720
00721 kreg->Set("FillkNNUseTrack", "no");
00722 kreg->Set("FillkNNTrim", "yes");
00723 kreg->Set("FillkNNTrimDelta", int(1000));
00724
00725 kreg->Set("FillkNNFilePath",baseConf+"/NuMuBar/data/knn.physics.notrack.far.daikon_04.dogwood1.L010z185i.root");
00726 if(fDetector == Detector::kNear)
00727 kreg->Set("FillkNNFilePath", baseConf+"/NuMuBar/data/knn.physics.notrack.near.daikon_04.dogwood1.L010z185i.root");
00728
00729
00730 kreg->Set("FillkNNKeyTruth", 19000);
00731
00732 kreg->Set("FillkNNKeyBase", 19400);
00733
00734 kreg->Set("FillkNNKeyList", "19008 19007 19004");
00735 kreg->Set("FillShortEventKeyPass",14001);
00736 kreg->Set("FillShortEventKeyPass2",14001);
00737 kreg->Set("FillkNN", "no");
00738
00739
00740 regJM->Set("EventkNN",*kreg);
00741
00742 regJM->Set("InterfaceConfigPath", baseConf+"/PhysicsNtuple/Config/ShortConfig2009.txt");
00743 regJM->Set("FillkNNFilePath", baseConf+"/NuMuBar/data/knn.physics.Alt.far.daikon_04.dogwood1.L010z185i.root");
00744 if(fDetector == Detector::kNear)
00745 regJM->Set("FillkNNFilePath", baseConf+"/NuMuBar/data/knn.physics.Alt.near.daikon_04.dogwood1.L010z185i.root");
00746
00747 fAnpInterfaceRO = new Anp::Interface();
00748 fAnpInterfaceRO->Config(*regRO);
00749
00750 fAnpInterfaceJM = new Anp::Interface();
00751 fAnpInterfaceJM->Config(*regJM);
00752 fkNNSet = true;
00753 }
00754
00755 fAnpInterfaceRO->FillSnarl(ntpManipulator->GetNtpStRecord());
00756 fAnpInterfaceJM->FillSnarl(ntpManipulator->GetNtpStRecord());
00757
00758 return;
00759
00760 }
00761
00762
00763
00764 void ANtpInfoObjectFillerNC::FillTrackInformation(ANtpRecoNtpManipulator *ntpManipulator,
00765 NtpSRTrack *ntpTrack,
00766 NtpSREvent *ntpEvent,
00767 ANtpTrackInfoNC *trackInfo,
00768 ANtpEventInfoNC *eventInfo)
00769 {
00770
00771 trackInfo->Reset();
00772 ANtpInfoObjectFiller::FillTrackInformation(ntpTrack, trackInfo);
00773
00774 if(trackInfo->totalStrips > 0){
00775 trackInfo->phPerStrip = trackInfo->pulseHeight/(1.*trackInfo->totalStrips);
00776 trackInfo->phPerPlane = trackInfo->pulseHeight/(1.*trackInfo->planes);
00777 }
00778
00779 NtpSRStrip *ntpStrip = 0;
00780
00781 Int_t plane = 0, strip = 0, stpCtr = 0;
00782 Float_t xPos = 0., yPos = 0., zPos = 0.;
00783 Int_t uXTalkStripCtr = 0;
00784 Int_t vXTalkStripCtr = 0;
00785
00786 Int_t numDigits = 0, numDoubleEndedStrips = 0;
00787 Int_t stripsInPlane[500], trkStripsInPlane[500];
00788
00789
00790
00791
00792 Double_t sumZT(0.); Double_t sumZ2(0.);
00793 Double_t sumT(0.); Double_t sumZ(0.);
00794 Int_t nHits(0);
00795
00796
00797
00798 ntpTrack = ntpManipulator->GetEventManipulator()->GetPrimaryTrackNS();
00799
00800 trackInfo->aTrkpass_ns =ntpTrack->fit.pass;
00801 trackInfo->aTrkph_ns =ntpTrack->ph.sigcor;
00802 trackInfo->aTrklen_ns =ntpTrack->plane.ntrklike;
00803 if(ntpTrack->plane.ntrklike>0) trackInfo->aTrkphperpl_ns = ntpTrack->ph.sigcor/ntpTrack->plane.ntrklike;
00804 if(ntpEvent->ph.sigcor>0) trackInfo->aTrkphper_ns = ntpTrack->ph.sigcor/ntpEvent->ph.sigcor;
00805 trackInfo->aTrkplu_ns =ntpTrack->plane.nu;
00806 trackInfo->aTrkplv_ns =ntpTrack->plane.nv;
00807 trackInfo->aTrkstp_ns =ntpTrack->nstrip;
00808 trackInfo->aTrkvtx_ns =ntpTrack->vtx.z;
00809
00810 for(Int_t ii = 0; ii < 500; ++ii){
00811 stripsInPlane[ii] = 0;
00812 trkStripsInPlane[ii] = 0;
00813 }
00814
00815
00816
00817 Float_t totalPH = 0.;
00818
00819 for(Int_t i = 0; i < ntpEvent->nstrip; ++i){
00820
00821 if(ntpEvent->stp[i] >= 0)
00822 ntpStrip = dynamic_cast<NtpSRStrip *>
00823 (fStripArray->At(ntpEvent->stp[i]));
00824 else continue;
00825 plane = (int)ntpStrip->plane;
00826 strip = ntpStrip->strip;
00827 totalPH = 0;
00828
00829 if( fDetector == Detector::kFar
00830 && !fStripIsXTalkEast[plane][strip] ) totalPH += ntpStrip->ph0.sigcor;
00831 if( !fStripIsXTalkWest[plane][strip] ) totalPH += ntpStrip->ph1.sigcor;
00832
00833 if(totalPH >= 90.){
00834 MSG("ANtpInfoObjectFillerNC", Msg::kDebug)
00835 << "plane " << plane << " strip " << ntpStrip->strip << " is ok "
00836 << fStripIsXTalkEast[plane][strip] << "/"
00837 << fStripIsXTalkWest[plane][strip] << " "
00838 << ntpStrip->ph0.sigcor << "/" << ntpStrip->ph1.sigcor << endl;
00839 ++stripsInPlane[(int)ntpStrip->plane];
00840 }
00841
00842 }
00843
00844
00845
00846
00847
00848
00849 for(Int_t i = 0; i < ntpTrack->nstrip; ++i){
00850
00851
00852 if (ntpTrack->stp[i] >= 0)
00853 ntpStrip = dynamic_cast<NtpSRStrip *>
00854 (fStripArray->At(ntpTrack->stp[i]));
00855 else continue;
00856 plane = (int)ntpStrip->plane;
00857 totalPH = 0;
00858
00859 if( fDetector > 0 && !fStripIsXTalkEast[plane][strip] ) totalPH += ntpStrip->ph0.sigcor;
00860 if( !fStripIsXTalkWest[plane][strip] ) totalPH += ntpStrip->ph1.sigcor;
00861
00862 if(totalPH >= 90.) ++trkStripsInPlane[(int)ntpStrip->plane];
00863
00864 if(fDetector == Detector::kFar
00865 && fStripIsXTalkEast[plane][strip] && fStripIsXTalkWest[plane][strip]){
00866 if((int)ntpStrip->planeview == kU) ++uXTalkStripCtr;
00867 else if((int)ntpStrip->planeview == kV) ++vXTalkStripCtr;
00868 }
00869 else if(fDetector == Detector::kNear && fStripIsXTalkWest[plane][strip]){
00870 if((int)ntpStrip->planeview == kU) ++uXTalkStripCtr;
00871 else if((int)ntpStrip->planeview == kV) ++vXTalkStripCtr;
00872 }
00873
00874 if (ntpStrip->ph0.pe > 2.) {
00875 sumZ += ntpStrip->z;
00876 sumZ2 += ntpStrip->z * ntpStrip->z;
00877 sumZT += ntpStrip->z * ntpStrip->time0;
00878 sumT += ntpStrip->time0;
00879 ++nHits;
00880 }
00881
00882 if (ntpStrip->ph1.pe > 2.) {
00883 sumZ += ntpStrip->z;
00884 sumZ2 += ntpStrip->z * ntpStrip->z;
00885 sumZT += ntpStrip->z * ntpStrip->time1;
00886 sumT += ntpStrip->time1;
00887 ++nHits;
00888 }
00889
00890
00891 }
00892
00893 trackInfo->xTalkStrips = uXTalkStripCtr+vXTalkStripCtr;
00894
00895
00896
00897 trackInfo->dtdz = 0;
00898 if (nHits > 1){
00899 Double_t CovZT = sumZT - (sumZ*sumT)/(Float_t)nHits;
00900 Double_t VarZ = sumZ2 - (sumZ*sumZ)/(Float_t)nHits;
00901 if (VarZ > 0) trackInfo->dtdz = Munits::c_light * CovZT/VarZ;
00902 }
00903
00904
00905
00906
00907
00908 trackInfo->trackLikePlanes = 0;
00909 for(Int_t i = 0; i < 500; ++i){
00910
00911 if(stripsInPlane[i]>0 && trkStripsInPlane[i]>0
00912 && stripsInPlane[i]*0.9 <= trkStripsInPlane[i])
00913 ++trackInfo->trackLikePlanes;
00914
00915 }
00916
00917
00918 Int_t index = 0;
00919 Int_t nVPlanes = 0, nUPlanes = 0;
00920
00921 nVPlanes = ntpTrack->plane.nv;
00922 nUPlanes = ntpTrack->plane.nu;
00923
00924 trackInfo->uvAsymmetry = 0.;
00925 if(nVPlanes+nUPlanes>0)
00926 trackInfo->uvAsymmetry = TMath::Abs((1.*nUPlanes-1.*nVPlanes)/(1.*nVPlanes+1.*nUPlanes));
00927
00928 trackInfo->signalUseFraction = 0.;
00929 if(eventInfo->pulseHeight>0.)
00930 trackInfo->signalUseFraction = trackInfo->pulseHeight/eventInfo->pulseHeight;
00931
00932 trackInfo->planeUseFraction = 0.;
00933 if(eventInfo->planes>0)
00934 trackInfo->planeUseFraction = 1.*trackInfo->trackLikePlanes/(1.*eventInfo->planes);
00935
00936 trackInfo->totalStrips = 0;
00937
00938
00939 if(trackInfo->length<50.){
00940
00941
00942 trackInfo->totalStrips = (int)ntpTrack->nstrip;
00943 numDoubleEndedStrips = 0;
00944
00945 stpCtr = 0;
00946
00947 for(Int_t ns = 0; ns < trackInfo->totalStrips; ++ns){
00948
00949
00950 if (ntpTrack->stp[ns] >= 0) index = ntpTrack->stp[ns];
00951 else continue;
00952
00953 ntpStrip = dynamic_cast<NtpSRStrip *>(fStripArray->At(index));
00954 plane = (int)ntpStrip->plane;
00955 strip = (int)ntpStrip->strip;
00956 numDigits = (int)ntpStrip->ndigit;
00957
00958 xPos = ntpTrack->stpx[ns];
00959 yPos = ntpTrack->stpy[ns];
00960 zPos = ntpTrack->stpz[ns];
00961
00962 if(numDigits == 2){
00963
00964 ++numDoubleEndedStrips;
00965
00966 }
00967 }
00968
00969 trackInfo->twoEndStripFraction = 0.;
00970 if(trackInfo->totalStrips>0)
00971 trackInfo->twoEndStripFraction = (1.*numDoubleEndedStrips)/(1.*trackInfo->totalStrips);
00972
00973 }
00974
00975
00976 trackInfo->kNN = fAnpInterfaceRO->GetVar("knn_pid",ntpEvent);
00977 trackInfo->numScintPlanes = fAnpInterfaceRO->GetVar("knn_01",ntpEvent);
00978 trackInfo->meanSigCor = fAnpInterfaceRO->GetVar("knn_10",ntpEvent);
00979 trackInfo->meanLowStripDivHighStrip = fAnpInterfaceRO->GetVar("knn_20",ntpEvent);
00980 trackInfo->trackSigCorFraction = fAnpInterfaceRO->GetVar("knn_40",ntpEvent);
00981
00982 trackInfo->kNNShort = fAnpInterfaceJM->GetVar("shortid",ntpEvent);
00983 trackInfo->kNNEvent = fAnpInterfaceJM->GetVar("evtid",ntpEvent);
00984
00985 return;
00986 }
00987
00988
00989
00990 void ANtpInfoObjectFillerNC::FillShowerInformation(NtpSRShower *ntpShower,
00991 ANtpShowerInfoNC *showerInfo,
00992 ANtpEventInfoNC *eventInfo)
00993 {
00994
00995 showerInfo->Reset();
00996 ANtpInfoObjectFiller::FillShowerInformation(ntpShower, showerInfo);
00997
00998 if(showerInfo->totalStrips > 0){
00999 showerInfo->phPerStrip = showerInfo->pulseHeight/(1.*showerInfo->totalStrips);
01000 showerInfo->phPerPlane = showerInfo->pulseHeight/(1.*showerInfo->planes);
01001 }
01002
01003 NtpSRStrip *ntpStrip = 0;
01004
01005 showerInfo->energyGeV = ntpShower->ph.gev;
01006
01007 Int_t index = 0;
01008 Int_t numDigits = 0;
01009 Int_t numDoubleEndedStrips = 0;
01010
01011 showerInfo->planeUseFraction = -1.;
01012 showerInfo->signalUseFraction = -1.;
01013 if(eventInfo->pulseHeight>0.)
01014 showerInfo->signalUseFraction = showerInfo->pulseHeight/(eventInfo->pulseHeight);
01015 if(eventInfo->planes>0)
01016 showerInfo->planeUseFraction = 1.*showerInfo->planes/(1.*eventInfo->planes);
01017
01018 Int_t plane = 0;
01019 Int_t strip = 0;
01020
01021 showerInfo->xTalkStrips = 0;
01022 showerInfo->twoEndStripFraction = -1.;
01023
01024 Float_t sumTposU = 0.;
01025 Float_t sum2TposU = 0.;
01026 Float_t sumPHU = 0.;
01027 Float_t maxTposU = -10.;
01028 Float_t minTposU = 10.;
01029
01030 Float_t sumTposV = 0.;
01031 Float_t sum2TposV = 0.;
01032 Float_t sumPHV = 0.;
01033 Float_t maxTposV = -10.;
01034 Float_t minTposV = 10.;
01035
01036 Float_t sumPHsigcor(0.);
01037 Int_t nstrip = ntpShower->nstrip;
01038
01039 MSG("ANtpInfoObjectFillerNC", Msg::kDebug) << "shower strips = "
01040 << ntpShower->nstrip << endl;
01041
01042 showerInfo->aShwdig_ns = ntpShower->ndigit;
01043 showerInfo->aShwplu_ns = ntpShower->plane.nu;
01044 showerInfo->aShwplv_ns = ntpShower->plane.nv;
01045
01046
01047
01048
01049
01050
01051
01052 for(Int_t ns = 0; ns < ntpShower->nstrip; ++ns){
01053
01054
01055 index = ntpShower->stp[ns];
01056
01057
01058 ntpStrip = dynamic_cast<NtpSRStrip *>(fStripArray->At(index));
01059 plane = (int)ntpStrip->plane;
01060 strip = ntpStrip->strip;
01061
01062
01063
01064
01065 sumPHsigcor += (ntpStrip->ph0.sigcor + ntpStrip->ph1.sigcor);
01066
01067
01068 numDigits = (int)ntpStrip->ndigit;
01069 if(numDigits == 2) ++numDoubleEndedStrips;
01070
01071 if(fDetector == Detector::kNear){
01072 if(fStripIsXTalkWest[plane][strip]) ++showerInfo->xTalkStrips;
01073 }
01074 else if(fDetector == Detector::kFar){
01075 if(fStripIsXTalkEast[plane][strip] && fStripIsXTalkWest[plane][strip])
01076 ++showerInfo->xTalkStrips;
01077 }
01078
01079 if(ntpStrip->ndigit == 1) continue;
01080 Float_t stpPH(ntpStrip->ph0.sigcor + ntpStrip->ph1.sigcor);
01081 if (ntpStrip->planeview == PlaneView::kU){
01082 maxTposU = TMath::Max(maxTposU, ntpStrip->tpos);
01083 minTposU = TMath::Min(minTposU, ntpStrip->tpos);
01084 sumTposU += (ntpStrip->tpos)*stpPH;
01085 sum2TposU += (ntpStrip->tpos * ntpStrip->tpos) * stpPH;
01086 sumPHU += stpPH;
01087 }
01088 else if (ntpStrip->planeview == PlaneView::kV){
01089 maxTposV = TMath::Max(maxTposV, ntpStrip->tpos);
01090 minTposV = TMath::Min(minTposV, ntpStrip->tpos);
01091 sumTposV += ntpStrip->tpos * stpPH;
01092 sum2TposV += (ntpStrip->tpos * ntpStrip->tpos) * stpPH;
01093 sumPHV +=stpPH;
01094 }
01095
01096 }
01097
01098 MSG("ANtpInfoObjectFillerNC", Msg::kDebug)
01099 << "total strips = " << showerInfo->totalStrips << endl;
01100
01101 if(showerInfo->totalStrips>0)
01102 showerInfo->twoEndStripFraction = (1.*numDoubleEndedStrips)/(1.*showerInfo->totalStrips);
01103
01104 if(sumPHU>0){
01105 Float_t centroidU = sumTposU;
01106 if(sumPHU > 0.) centroidU /= sumPHU;
01107 else centroidU = 0.;
01108 showerInfo->transverseRMSU = TMath::Sqrt(TMath::Abs(sum2TposU/sumPHU - centroidU*centroidU));
01109 }
01110 else showerInfo->transverseRMSU = 0.;
01111
01112 if(sumPHV>0){
01113 Float_t centroidV = sumTposV;
01114 if(sumPHV > 0.) centroidV /= sumPHV;
01115 else centroidV = 0.;
01116 showerInfo->transverseRMSV = TMath::Sqrt(TMath::Abs(sum2TposV/sumPHV - centroidV*centroidV));
01117 }
01118 else showerInfo->transverseRMSV = 0.;
01119
01120
01121
01122
01123
01124 Float_t meanPHsigcor = sumPHsigcor/nstrip;
01125 Float_t secondMoment(0.);
01126 Float_t fourthMoment(0.);
01127
01128 for (Int_t stp_index_shw = 0; stp_index_shw < nstrip ; ++stp_index_shw){
01129 Int_t stp_index = ntpShower->stp[stp_index_shw];
01130 NtpSRStrip* strip =
01131 dynamic_cast<NtpSRStrip *>(fStripArray->At(stp_index));
01132 Float_t stpPH(strip->ph0.sigcor + strip->ph1.sigcor);
01133 secondMoment += TMath::Power((stpPH - meanPHsigcor),2);
01134 fourthMoment += TMath::Power((stpPH - meanPHsigcor),4);
01135 }
01136 Float_t sigma = secondMoment/(nstrip-1);
01137
01138
01139 showerInfo->phKurtosis = fourthMoment/(nstrip*sigma*sigma)-3.0;
01140
01141 Int_t emLikeStrips(0);
01142 Bool_t cluFailure(false);
01143 for (Int_t clu_index_shw = 0; clu_index_shw < ntpShower->ncluster ;
01144 ++clu_index_shw){
01145 Int_t clu_index = ntpShower->clu[clu_index_shw];
01146
01147 if (clu_index<0) {
01148 MSG("ANtpInfoObjectFillerNC",Msg::kError)
01149 << "Array index for cluster less than 0. This should never happen!"
01150 << endl;
01151 cluFailure = true;
01152 continue;
01153 }
01154
01155 NtpSRCluster* cluster =
01156 dynamic_cast<NtpSRCluster *>(fClusterArray->At(clu_index));
01157 if (cluster->id == 0)
01158 emLikeStrips += cluster->nstrip;
01159 }
01160 if (!cluFailure)
01161 showerInfo->emFrac = (1.0* emLikeStrips)/ntpShower->nstrip;
01162 else
01163 showerInfo->emFrac = -1;
01164
01165
01166
01167 return;
01168 }
01169
01170
01171
01172 void ANtpInfoObjectFillerNC::FillCrossOverInformation(NtpSRTrack *ntpTrack,
01173 NtpSRShower *ntpShower,
01174 NtpSREvent *ntpEvent,
01175 ANtpRecoNtpManipulator *ntpManipulator,
01176 ANtpTrackInfoNC *trackInfo,
01177 ANtpShowerInfoNC *showerInfo,
01178 ANtpEventInfoNC *eventInfo)
01179 {
01180
01181 Int_t indexShowerStrip = 0;
01182 Int_t indexTrackStrip = 0;
01183
01184 eventInfo->trackStripsInShower = 0;
01185 eventInfo->trackExtension = 0;
01186
01187 if(ntpShower && ntpTrack && trackInfo->length < 50.){
01188 for(Int_t ns = 0; ns < ntpShower->nstrip; ++ns){
01189
01190
01191 indexShowerStrip = ntpShower->stp[ns];
01192
01193
01194 for(Int_t nt = 0; nt < ntpTrack->nstrip; ++nt){
01195
01196 indexTrackStrip = ntpTrack->stp[nt];
01197
01198
01199
01200 if(indexTrackStrip == indexShowerStrip) ++eventInfo->trackStripsInShower;
01201
01202 }
01203
01204 }
01205
01206 if(ANtpDefaultValue::IsDefault(showerInfo->planes) ||
01207 ANtpDefaultValue::IsDefault(trackInfo->planes) )
01208 eventInfo->trackExtension = ANtpDefaultValue::kInt;
01209 else
01210 eventInfo->trackExtension = showerInfo->planes - trackInfo->planes;
01211 }
01212
01213
01214 Int_t ssind,ssplind;
01215 Double_t ssphtot;
01216 Bool_t foundshw,foundtrk;
01217
01218 Int_t planes = ntpEvent->plane.beg;
01219
01220 Double_t ph3,ph6,phlast,phcommon,phnowere,hitcommon,hitnowere;
01221 ph3=0.;ph6=0.;phlast=0.;phcommon=0.;phnowere=0.;hitcommon=0.;hitnowere=0.;
01222
01223 NtpSRTrack *ntpTrackNS = ntpManipulator->GetEventManipulator()->GetPrimaryTrackNS();
01224
01225
01226 for(Int_t evss=0;evss<ntpEvent->nstrip;evss++){
01227 Int_t stp_index=((ntpEvent->stp)[evss]);
01228
01229 if(stp_index!=-1){
01230
01231 NtpSRStrip *ntpStrip = dynamic_cast<NtpSRStrip *>(fStripArray->At(stp_index));
01232 if(!ntpStrip) continue;
01233 ssind = ntpStrip->strip;
01234 ssplind = ntpStrip->plane;
01235 ssphtot = ntpStrip->ph1.sigcor+ntpStrip->ph0.sigcor;
01236
01237 foundshw=false;
01238 foundtrk=false;
01239
01240 Double_t phstrips=0;
01241 Double_t phstript=0;
01242
01243 Int_t sshwind,sshwplind;
01244 Double_t sshwphtot;
01245
01246 if(ntpEvent->nshower > 0) {
01247 for(Int_t jj=0;jj<ntpShower->nstrip;jj++){
01248 Int_t stp_indexs=((ntpShower->stp)[jj]);
01249
01250 if(stp_indexs!=-1){
01251
01252 NtpSRStrip *ntpStrip1 = dynamic_cast<NtpSRStrip *>(fStripArray->At(stp_indexs));
01253 if(!foundshw){
01254 sshwind =ntpStrip1->strip;
01255 sshwplind =ntpStrip1->plane;
01256 sshwphtot =ntpStrip1->ph1.sigcor+ntpStrip1->ph0.sigcor;
01257 if(sshwind==ssind && sshwplind==ssplind) {
01258 foundshw=true;
01259 phstrips=sshwphtot;
01260 }
01261 }
01262 }
01263 }
01264 }
01265
01266 Int_t strkind,strkplind;
01267 Double_t strkphtot;
01268
01269 if(ntpEvent->ntrack > 0) {
01270
01271 for(Int_t jj=0;jj<ntpTrackNS->nstrip;jj++){
01272 Int_t stp_indext=((ntpTrackNS->stp)[jj]);
01273
01274 if(stp_indext!=-1){
01275
01276 NtpSRStrip *ntpStript = dynamic_cast<NtpSRStrip *>(fStripArray->At(stp_indext));
01277
01278 if(!foundtrk){
01279 strkind =ntpStript->strip;
01280 strkplind=ntpStript->plane;
01281 strkphtot=ntpStript->ph1.sigcor+ntpStript->ph0.sigcor;
01282 if(strkind==ssind && strkplind==ssplind) {
01283 foundtrk=true;
01284 phstript=strkphtot;
01285 }
01286 }
01287
01288 }
01289 }
01290 }
01291
01292 if(foundshw && foundtrk) {
01293 hitcommon= hitcommon+1;
01294 phcommon = phcommon+phstrips+phstript;
01295 }
01296 if(!foundshw && ! foundtrk) {
01297 hitnowere=hitnowere+1;
01298 phnowere=phnowere+ssphtot;
01299 }
01300 if(ssplind>=planes && ssplind<=(planes+3)){
01301 ph3=ph3+ssphtot;
01302 }
01303 else if(ssplind>(planes+3) && ssplind<=(planes+6)){
01304 ph6=ph6+ssphtot;
01305 }
01306 else {
01307 phlast=phlast+ssphtot;
01308 }
01309
01310 }
01311
01312 }
01313
01314
01315 eventInfo->aPh3_ns = ph3;
01316 eventInfo->aPh6_ns = ph6;
01317 eventInfo->aPhlast_ns = phlast;
01318 eventInfo->aPhcommon_ns = phcommon;
01319
01320
01321 return;
01322 }
01323
01324
01325
01326 void ANtpInfoObjectFillerNC::FillEventInformation(ANtpRecoNtpManipulator *ntpManipulator,
01327 NtpSREvent *ntpEvent,
01328 ANtpEventInfoNC *eventInfo)
01329 {
01330
01331 eventInfo->Reset();
01332 ANtpInfoObjectFiller::FillEventInformation(ntpManipulator,
01333 ntpEvent,
01334 eventInfo);
01335
01336 NtpSRStrip *ntpStrip = 0;
01337
01338 eventInfo->lengthInPlanes = eventInfo->endPlane - eventInfo->begPlane + 1;
01339 eventInfo->lateBucketPHFraction = ntpEvent->bleach.lateBucketPHFraction;
01340 eventInfo->timeWeightedPHFraction = ntpEvent->bleach.timeWeightedPHFraction;
01341 eventInfo->straightPHFraction = ntpEvent->bleach.straightPHFraction;
01342 eventInfo->fixedWindowPH = ntpEvent->bleach.fixedWindowPH;
01343
01344
01345 Float_t planePH[500] = {0.};
01346 eventInfo->totalStrips = ntpEvent->nstrip;
01347 eventInfo->passStrips = 0;
01348 eventInfo->modifiedPH = 0.;
01349 if(eventInfo->totalStrips > 0){
01350 eventInfo->phPerStrip = eventInfo->pulseHeight/(1.*eventInfo->totalStrips);
01351 eventInfo->phPerPlane = eventInfo->pulseHeight/(1.*eventInfo->planes);
01352 }
01353
01354 double maxStripTime = -1.e12;
01355 double minStripTime = 1.e12;
01356
01357
01358
01359
01360
01361
01362 const UInt_t nPlanes(eventInfo->planes);
01363 vector<Float_t> planeAdc(nPlanes+12,0);
01364
01365 const vector<Float_t>::const_iterator pAdc_begin(planeAdc.begin());
01366 const vector<Float_t>::const_iterator pAdc_last(pAdc_begin + nPlanes);
01367
01368
01369 Double_t t_1stStp(10.);
01370 Double_t t_lastStp(0.);
01371 Double_t t_sumStp(0.);
01372 Double_t t_sum2Stp(0.);
01373 Int_t nStpHit( ntpEvent->nstrip );
01374
01375
01376 for(Int_t ns = 0; ns < ntpEvent->nstrip; ++ns){
01377
01378
01379 if (ntpEvent->stp[ns] >= 0)
01380 ntpStrip = dynamic_cast<NtpSRStrip *>
01381 (fStripArray->At(ntpEvent->stp[ns]));
01382 else continue;
01383 if(ntpStrip->plane>0 && ntpStrip->plane<486){
01384 planePH[ntpStrip->plane] += ntpStrip->ph0.sigcor + ntpStrip->ph1.sigcor;
01385
01386 if( fDetector == Detector::kFar
01387 && (!fStripIsXTalkWest[ntpStrip->plane][ntpStrip->strip]
01388 || !fStripIsXTalkEast[ntpStrip->plane][ntpStrip->strip]))
01389 ++eventInfo->passStrips;
01390 else if( fDetector == Detector::kNear
01391 && !fStripIsXTalkWest[ntpStrip->plane][ntpStrip->strip])
01392 ++eventInfo->passStrips;
01393
01394 if(ntpStrip->ph1.raw > 200.) eventInfo->modifiedPH += ntpStrip->ph1.raw;
01395
01396 if(ntpStrip->time1 > maxStripTime) maxStripTime = ntpStrip->time1;
01397 if(ntpStrip->time1 < minStripTime) minStripTime = ntpStrip->time1;
01398
01399 }
01400
01401
01402
01403
01404 Int_t pln_index = ntpStrip->plane - eventInfo->begPlane;
01405 if (pln_index < (Int_t) nPlanes && pln_index >= 0) {
01406 planeAdc[pln_index] += ntpStrip->ph0.sigcor;
01407 planeAdc[pln_index] += ntpStrip->ph1.sigcor;
01408 }
01409
01410
01411 Double_t t_stpHitA = TMath::Max(ntpStrip->time0 , ntpStrip->time1);
01412 Double_t t_stpHitB = TMath::Min(ntpStrip->time0 , ntpStrip->time1);
01413
01414 t_1stStp = TMath::Min(t_1stStp, t_stpHitA);
01415 t_lastStp = TMath::Max(t_lastStp, t_stpHitA);
01416 t_sumStp += t_stpHitA;
01417 t_sum2Stp += t_stpHitA * t_stpHitA;
01418
01419 if (t_stpHitB > 0 ){
01420 t_1stStp = TMath::Min(t_1stStp, t_stpHitB);
01421 t_lastStp = TMath::Max(t_lastStp, t_stpHitB);
01422 t_sumStp += t_stpHitB;
01423 t_sum2Stp += t_stpHitB * t_stpHitB;
01424 ++nStpHit;
01425 }
01426
01427
01428 }
01429
01430
01431
01432 eventInfo->stripTimeMean = t_sumStp / nStpHit;
01433 Float_t t_variance = t_sum2Stp / nStpHit
01434 - (t_sumStp / nStpHit) * (t_sumStp / nStpHit);
01435
01436 eventInfo->stripTimeRMS = TMath::Sqrt(TMath::Abs(t_variance));
01437
01438 eventInfo->stripTime1st = t_1stStp;
01439 eventInfo->stripTimelast = t_lastStp;
01440
01441
01442 eventInfo->eventDuration = maxStripTime - minStripTime;
01443
01444 if(fDetector == Detector::kNear) FindEarlyActivityWeight(ntpEvent, eventInfo);
01445
01446
01447
01448
01449
01450 TH1F planePHDist("planePHDist", "", 400, 0., 400.);
01451 eventInfo->maxPlanePH = 0.;
01452 eventInfo->maxPHIn3Planes = 0.;
01453 eventInfo->maxPHIn6Planes = 0.;
01454 eventInfo->maxPHIn9Planes = 0.;
01455 eventInfo->maxPHIn12Planes = 0.;
01456
01457 for(Int_t i = eventInfo->begPlane; i <= eventInfo->endPlane; ++i){
01458 if(planePH[i] > eventInfo->maxPlanePH) eventInfo->maxPlanePH = planePH[i];
01459 if(planePH[i] > 0.) planePHDist.Fill(planePH[i]*1.e-4);
01460 }
01461
01462 Float_t phInWindow = 0.;
01463 for(Int_t windowSize = 3; windowSize<14; windowSize +=3){
01464 for(Int_t i = eventInfo->begPlane; i<=eventInfo->endPlane; ++i){
01465 phInWindow = 0.;
01466 for(Int_t j = 0; j<windowSize && j+i<=eventInfo->endPlane; ++j){
01467 phInWindow += planePH[i+j];
01468 }
01469 if(windowSize == 3 && eventInfo->maxPHIn3Planes<phInWindow)
01470 eventInfo->maxPHIn3Planes = phInWindow;
01471 else if(windowSize == 6 && eventInfo->maxPHIn6Planes<phInWindow)
01472 eventInfo->maxPHIn6Planes = phInWindow;
01473 else if(windowSize == 9 && eventInfo->maxPHIn9Planes<phInWindow)
01474 eventInfo->maxPHIn9Planes = phInWindow;
01475 else if(windowSize == 12 && eventInfo->maxPHIn12Planes<phInWindow)
01476 eventInfo->maxPHIn12Planes = phInWindow;
01477 }
01478
01479 }
01480
01481 eventInfo->pulseHeightRms = planePHDist.GetRMS();
01482
01483
01484
01485
01486
01487 eventInfo->triPlane1PH = accumulate(pAdc_begin, pAdc_begin+3, 0.0 );
01488 eventInfo->triPlane2PH = accumulate(pAdc_begin+3, pAdc_begin+6, 0.0 );
01489 if (nPlanes>6)
01490 eventInfo->triPlaneOverPH = accumulate(pAdc_begin+6, pAdc_last, 0.0 );
01491 else
01492 eventInfo->triPlaneOverPH = 0.0;
01493
01494
01495 eventInfo->eventSummaryPlanes = ntpManipulator->GetSnarlEventSummary().plane.n;
01496
01497
01498
01499
01500 return;
01501 }
01502
01503
01504
01505 void ANtpInfoObjectFillerNC::FillMCInformation(NtpMCTruth *ntpMCTruth,
01506 TClonesArray *stdHepArray,
01507 ANtpTruthInfoBeam *truthInfo)
01508 {
01509 truthInfo->Reset();
01510 ANtpInfoObjectFiller::FillMCTruthInformation(ntpMCTruth, truthInfo);
01511 ANtpInfoObjectFillerBeam::FillBeamMCTruthInformation(ntpMCTruth, stdHepArray, truthInfo);
01512
01513 return;
01514 }
01515
01516
01517 void ANtpInfoObjectFillerNC::FillPlanePixelSignalArrays(NtpSREvent *ntpEvent)
01518 {
01519 MSG("ANtpInfoObjectFillerNC", Msg::kDebug) << "in FillPlanePixelSignalArrays" << endl;
01520
01521 NtpSRStrip *ntpStrip = 0;
01522
01523 Int_t pixel = 0;
01524 Int_t plane = 0;
01525
01526
01527
01528
01529
01530
01531 for(Int_t ns = 0; ns < ntpEvent->nstrip; ++ns){
01532
01533
01534 if (ntpEvent->stp[ns] >= 0)
01535 ntpStrip = dynamic_cast<NtpSRStrip *>
01536 (fStripArray->At(ntpEvent->stp[ns]));
01537 else continue;
01538 plane = ntpStrip->plane;
01539
01540 MSG("ANtpInfoObjectFillerNC", Msg::kDebug)
01541 << "plane " << plane << " strip " << ntpStrip->strip
01542 << " pmt1 " << ntpStrip->pmtindex1 << endl;
01543 MSG("ANtpInfoObjectFillerNC", Msg::kDebug)
01544 << " pmt0 " << ntpStrip->pmtindex0 << endl;
01545
01546 fPlaneToPMTMapWest[plane] = ntpStrip->pmtindex1;
01547 if(fDetector == Detector::kFar) fPlaneToPMTMapEast[plane] = ntpStrip->pmtindex0;
01548
01549 }
01550
01551 MSG("ANtpInfoObjectFillerNC", Msg::kDebug)
01552 << "finish filling plane to pmt maps" << endl;
01553
01554
01555 for(Int_t ns = 0; ns < ntpEvent->nstrip; ++ns){
01556
01557
01558 if(ntpEvent->stp[ns] >= 0)
01559 ntpStrip = dynamic_cast<NtpSRStrip *>
01560 (fStripArray->At(ntpEvent->stp[ns]));
01561 else continue;
01562 plane = ntpStrip->plane;
01563
01564 MSG("ANtpInfoObjectFillerNC", Msg::kDebug)
01565 << "plane " << plane << " strip " << ntpStrip->strip
01566 << " detector " << fDetector << endl;
01567
01568
01569 if(fDetector == Detector::kFar){
01570
01571
01572 pixel = fStripToPixelMapEast[ntpStrip->strip];
01573
01574 if(ntpStrip->pmtindex0 == fPlaneToPMTMapEast[plane])
01575 fPlanePixelEastSignal[plane][pixel][0] = ntpStrip->ph0.sigcor;
01576 else
01577 fPlanePixelEastSignal[plane][pixel][1] = ntpStrip->ph0.sigcor;
01578
01579
01580 pixel = fStripToPixelMapWest[ntpStrip->strip];
01581 if(ntpStrip->pmtindex1 == fPlaneToPMTMapWest[plane])
01582 fPlanePixelWestSignal[plane][pixel][0] = ntpStrip->ph1.sigcor;
01583 else
01584 fPlanePixelWestSignal[plane][pixel][1] = ntpStrip->ph1.sigcor;
01585
01586 }
01587 else if(fDetector == Detector::kNear){
01588
01589
01590
01591
01592 MSG("ANtpInfoObjectFillerNC", Msg::kDebug)
01593 << "in near detector algorithm "
01594 << " plane " << plane
01595 << " coverage = " << fPlaneCoverage[plane]
01596 << " view = " << (int)ntpStrip->planeview << endl;
01597
01598 if(fPlaneCoverage[plane] == kNearFull){
01599
01600 if((int)ntpStrip->planeview == kU){
01601 pixel = fStripToPixelMapNearU2[ntpStrip->strip];
01602 if(pixel < 0) pixel = fStripToPixelMapNearU3[ntpStrip->strip];
01603 if(pixel < 0)
01604 MSG("ANtpInfoObjectFillerNC", Msg::kDebug)
01605 << "cannot find pixel for strip " << ntpStrip->strip
01606 << " plane " << plane << endl;
01607
01608 }
01609 else if((int)ntpStrip->planeview == kV){
01610 pixel = fStripToPixelMapNearV2[ntpStrip->strip];
01611 if(pixel < 0) pixel = fStripToPixelMapNearV3[ntpStrip->strip];
01612 if(pixel < 0)
01613 MSG("ANtpInfoObjectFillerNC", Msg::kDebug)
01614 << "cannot find pixel for strip " << ntpStrip->strip
01615 << " plane " << plane << endl;
01616
01617 }
01618
01619 }
01620 else if(fPlaneCoverage[plane] == kNearPartial){
01621
01622 if((int)ntpStrip->planeview == kU)
01623 pixel = fStripToPixelMapNearU1[ntpStrip->strip];
01624 else if((int)ntpStrip->planeview == kV)
01625 pixel = fStripToPixelMapNearV1[ntpStrip->strip];
01626
01627 }
01628
01629 if(ntpStrip->pmtindex1 == fPlaneToPMTMapWest[plane])
01630 fPlanePixelWestSignal[plane][pixel][0] = ntpStrip->ph1.sigcor;
01631 else
01632 fPlanePixelWestSignal[plane][pixel][1] = ntpStrip->ph1.sigcor;
01633
01634
01635 }
01636 }
01637
01638 return;
01639 }
01640
01641
01642 void ANtpInfoObjectFillerNC::FindEarlyActivityWeight(NtpSREvent *ntpEvent,
01643 ANtpEventInfoNC *eventInfo)
01644 {
01645
01646
01647 double earliestEventTime = 1.e20;
01648 map<int,int> eventPlanes;
01649 NtpSRStrip *strip = 0;
01650
01651 for(int i = 0; i < ntpEvent->nstrip; ++i){
01652 if (ntpEvent->stp[i] >= 0)
01653 strip = dynamic_cast<NtpSRStrip *>
01654 (fStripArray->At(ntpEvent->stp[i]));
01655 else continue;
01656 if(strip->time1 < earliestEventTime) earliestEventTime = strip->time1;
01657 if(eventPlanes.find(strip->pmtindex1) == eventPlanes.end()) eventPlanes[strip->pmtindex1] = 1;
01658
01659 }
01660
01661
01662
01663
01664
01665
01666
01667
01668
01669
01670
01671
01672
01673
01674 double weightSum = 0.;
01675 for(int i = 0; i < fStripArray->GetLast()+1; ++i){
01676
01677 strip = dynamic_cast<NtpSRStrip *>(fStripArray->At(i));
01678
01679
01680 if(1.e9*(earliestEventTime - strip->time1) > 0.
01681 && 1.e9*(earliestEventTime - strip->time1) < 1000.*1.5
01682 && eventPlanes.find(strip->pmtindex1) != eventPlanes.end()){
01683 weightSum += strip->ph1.sigcor*TMath::Exp(-1.e9*(earliestEventTime - strip->time1)/700.);
01684
01685 MSG("ANtpInfoObjectFillerNC", Msg::kDebug)
01686 << strip->pmtindex1 << " " << strip->plane
01687 << " " << strip->ph1.sigcor
01688 << " " << 1.e9*(earliestEventTime-strip->time1)
01689 << " " << weightSum
01690 << " " << eventInfo->pulseHeight << endl;
01691
01692 }
01693 }
01694
01695 eventInfo->earlyWeightedADC = weightSum;
01696
01697 return;
01698 }
01699
01700
01701 Float_t ANtpInfoObjectFillerNC::FindNearestNeighborPixelSignal(Int_t plane, Int_t pixel, Int_t pmt,
01702 Float_t planePixelSignal[][64][2])
01703 {
01704 MSG("ANtpInfoObjectFillerNC", Msg::kDebug)
01705 << "in FindNearestNeighborPixelSignal" << endl;
01706
01707 Float_t sumNeighborSignal = 0.;
01708
01709 if(fDetector == Detector::kFar){
01710 if(pixel == 5 || pixel == 6 || pixel == 9 || pixel == 10)
01711 sumNeighborSignal = (planePixelSignal[plane][pixel-4][pmt] + planePixelSignal[plane][pixel+4][pmt]
01712 + planePixelSignal[plane][pixel-1][pmt] + planePixelSignal[plane][pixel+1][pmt]);
01713 else if(pixel < 3 && pixel > 0)
01714 sumNeighborSignal = (planePixelSignal[plane][pixel-1][pmt] + planePixelSignal[plane][pixel+1][pmt]
01715 + planePixelSignal[plane][pixel+4][pmt]);
01716 else if(pixel < 15 && pixel > 12)
01717 sumNeighborSignal = (planePixelSignal[plane][pixel-1][pmt] + planePixelSignal[plane][pixel+1][pmt]
01718 + planePixelSignal[plane][pixel-4][pmt]);
01719 else if(pixel > 3 && pixel < 15 && (pixel+1)%4 == 0)
01720 sumNeighborSignal = (planePixelSignal[plane][pixel-4][pmt] + planePixelSignal[plane][pixel+4][pmt]
01721 + planePixelSignal[plane][pixel-1][pmt]);
01722
01723 else if(pixel > 3 && pixel < 12 && pixel%4 == 0)
01724 sumNeighborSignal = (planePixelSignal[plane][pixel-4][pmt] + planePixelSignal[plane][pixel+4][pmt]
01725 + planePixelSignal[plane][pixel+1][pmt]);
01726 else if(pixel == 0)
01727 sumNeighborSignal = planePixelSignal[plane][pixel+4][pmt] + planePixelSignal[plane][pixel+1][pmt];
01728
01729 else if(pixel == 3)
01730 sumNeighborSignal = planePixelSignal[plane][pixel+4][pmt] + planePixelSignal[plane][pixel-1][pmt];
01731
01732 else if(pixel == 12)
01733 sumNeighborSignal = planePixelSignal[plane][pixel-4][pmt] + planePixelSignal[plane][pixel+1][pmt];
01734
01735 else if(pixel == 15)
01736 sumNeighborSignal = planePixelSignal[plane][pixel-4][pmt] + planePixelSignal[plane][pixel-1][pmt];
01737
01738 }
01739 else if(fDetector == Detector::kNear){
01740
01741 if( (pixel > 8 && pixel < 15) || (pixel > 16 && pixel < 23)
01742 || (pixel > 24 && pixel < 31) || (pixel > 32 && pixel < 39)
01743 || (pixel > 40 && pixel < 47) || (pixel > 48 && pixel < 55) )
01744 sumNeighborSignal = (planePixelSignal[plane][pixel-8][pmt] + planePixelSignal[plane][pixel+8][pmt]
01745 + planePixelSignal[plane][pixel-1][pmt] + planePixelSignal[plane][pixel+1][pmt]);
01746 else if( pixel > 0 && pixel < 7)
01747 sumNeighborSignal = (planePixelSignal[plane][pixel-1][pmt] + planePixelSignal[plane][pixel+1][pmt]
01748 + planePixelSignal[plane][pixel+8][pmt]);
01749 else if( pixel > 56 && pixel < 63)
01750 sumNeighborSignal = (planePixelSignal[plane][pixel-1][pmt] + planePixelSignal[plane][pixel+1][pmt]
01751 + planePixelSignal[plane][pixel-8][pmt]);
01752 else if( pixel > 8 && pixel < 63 && (pixel+1)%8==0 )
01753 sumNeighborSignal = (planePixelSignal[plane][pixel-8][pmt] + planePixelSignal[plane][pixel+8][pmt]
01754 + planePixelSignal[plane][pixel-1][pmt]);
01755 else if( pixel > 7 && pixel < 56 && (pixel)%8==0 )
01756 sumNeighborSignal = (planePixelSignal[plane][pixel-8][pmt] + planePixelSignal[plane][pixel+8][pmt]
01757 + planePixelSignal[plane][pixel+1][pmt]);
01758 else if(pixel==0)
01759 sumNeighborSignal = planePixelSignal[plane][pixel+8][pmt] + planePixelSignal[plane][pixel+1][pmt];
01760 else if(pixel==7)
01761 sumNeighborSignal = planePixelSignal[plane][pixel+8][pmt] + planePixelSignal[plane][pixel-1][pmt];
01762 else if(pixel==56)
01763 sumNeighborSignal = planePixelSignal[plane][pixel-8][pmt] + planePixelSignal[plane][pixel+1][pmt];
01764 else if(pixel==63)
01765 sumNeighborSignal = planePixelSignal[plane][pixel-8][pmt] + planePixelSignal[plane][pixel-1][pmt];
01766 }
01767
01768 return sumNeighborSignal;
01769 }
01770
01771
01772 void ANtpInfoObjectFillerNC::FindXTalkStrips(NtpSREvent *ntpEvent)
01773 {
01774 MSG("ANtpInfoObjectFillerNC", Msg::kDebug) << "in FindXTalkStrips" << endl;
01775
01776 NtpSRStrip *ntpStrip = 0;
01777
01778 Int_t plane = 0;
01779 Int_t pixelEast = 0;
01780 Int_t pmtEast = 0;
01781 Int_t pixelWest = 0;
01782 Int_t pmtWest = 0;
01783
01784 for(Int_t ns = 0; ns < ntpEvent->nstrip; ++ns){
01785
01786
01787 if (ntpEvent->stp[ns] >= 0)
01788 ntpStrip = dynamic_cast<NtpSRStrip *>
01789 (fStripArray->At(ntpEvent->stp[ns]));
01790 else continue;
01791 plane = (int)ntpStrip->plane;
01792
01793 if(fDetector == Detector::kFar){
01794 pixelEast = fStripToPixelMapEast[ntpStrip->strip];
01795 pmtEast = 0;
01796 if(fPlaneToPMTMapEast[plane] != ntpStrip->pmtindex0) pmtEast = 1;
01797
01798 if(ntpStrip->ph0.sigcor < 0.1*FindNearestNeighborPixelSignal(plane, pixelEast,
01799 pmtEast, fPlanePixelEastSignal)
01800 || ntpStrip->ph0.sigcor < 90.){
01801 fStripIsXTalkEast[plane][ntpStrip->strip] = 1;
01802 MSG("ANtpInfoObjectFillerNC", Msg::kDebug)
01803 << "east plane " << plane << " strip " << ntpStrip->strip
01804 << " is xtalk " << ntpStrip->ph0.sigcor << endl;
01805 }
01806
01807 }
01808
01809 pixelWest = fStripToPixelMapWest[ntpStrip->strip];
01810 pmtWest = 0;
01811 if(fPlaneToPMTMapWest[plane] != ntpStrip->pmtindex1) pmtWest = 1;
01812
01813 if(ntpStrip->ph1.sigcor < 0.1*FindNearestNeighborPixelSignal(plane, pixelWest,
01814 pmtWest, fPlanePixelWestSignal)
01815 || ntpStrip->ph1.sigcor < 90.){
01816
01817 fStripIsXTalkWest[plane][ntpStrip->strip] = 1;
01818 MSG("ANtpInfoObjectFillerNC", Msg::kDebug)
01819 << "west plane " << plane << " strip " << ntpStrip->strip
01820 << " is xtalk " << ntpStrip->ph1.sigcor << endl;
01821 }
01822
01823 }
01824 return;
01825 }
01826