loon: ANtpInfoObjectFillerNC.cxx Source File

00001 
00002 //$Id: ANtpInfoObjectFillerNC.cxx,v 1.50 2010/07/28 20:51:36 rhatcher Exp $
00003 //
00004 //ANtpInfoObjectFiller.cxx
00005 //
00006 //B. Rebel 12/2005
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 // #include "AnalysisNtuples/Module/VHS.h"
00038 // #include "AnalysisNtuples/Module/VHSevent.h"
00039 
00040 #include "Conventions/Munits.h"
00041 #include "Conventions/PlaneView.h"
00042 #include "Conventions/ReleaseType.h"
00043 
00044 #include "Registry/Registry.h"
00045 
00046 //ROOT includes
00047 #include "TFile.h"
00048 #include "TH1.h"
00049 #include "TSystem.h"
00050 
00051 //C++/STL includes
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   // ANtpInfoObjectFillerBeam(Detector::kUnknown),
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   // Fills the strip to pixel map. Call this when you pick a detector
00111 
00112   if(fDetector  == Detector::kFar){
00113 
00114     //Define the FD strip to pixel maps
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     // Fill the strip to pixe map data members
00134     for(Int_t i = 0; i < 192; ++i){
00135       fStripToPixelMapWest[i] = stripToPixelWest[i];
00136       fStripToPixelMapEast[i] = stripToPixelEast[i];
00137     }
00138 
00139     // set plane coverages
00140     for(Int_t i = 0; i < 486; ++i) fPlaneCoverage[i] = kComplete;
00141   } // end filling strip to pixel map for far detector
00142 
00143   else if(fDetector == Detector::kNear){
00144 
00145     // Define the pixel to strip maps
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     //now fill the strip to pixel map data members
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     }//end fill strip to pixel maps
00218 
00219     //fill the array of plane coverages
00220     for(Int_t i = 0; i < 282; ++i){
00221       if( i%5 == 1 ) fPlaneCoverage[i] = kNearFull;
00222       else fPlaneCoverage[i] = kNearPartial;
00223     }
00224   }//end filling strip to pixel map for near detector
00225 
00226   MSG("ANtpInfoObjectFillerNC", Msg::kDebug)
00227     << "ANtpInfoObjectFillerNC::Constructor" << endl;
00228 
00229 }
00230 
00231 //----------------------------------------------------------------------
00232 void ANtpInfoObjectFillerNC::SetDetector(Detector::Detector_t detector)
00233 {
00234   // Set the detector. Then fill the strip to pixel maps accordingly
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   //get event.  the call to ANtpEventManipulator::SetEventInSnarl sets the
00264   //NtpSREvent data member in that object to the current event so that
00265   //you can later call for the primary track, shower, etc.
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   //get the primary shower for the event - if no shower is present it
00274   //returns 0
00275   ntpShower = ntpManipulator->GetEventManipulator()->GetPrimaryShower();
00276   if(ntpShower) FillShowerInformation(ntpShower, showerInfo, eventInfo);
00277 
00278   //get the primary track for the event - if no track is present it
00279   //returns 0
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   // Get the truth information for this track if it is MC
00288   // Get best neu match from Truthhelper
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   // added by Tobi
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   //Correct cedar vertex bug
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){ // Take fixed vertex if possible
00320       vtxPlane = vtxf.VtxPlane();
00321       vtxZ = vtxf.VtxZ();
00322     }
00323   } // end if cedar
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     //fill longitudnal energy deposition histogram                                                                   
00347     if(deltaplanes>=0.)
00348       lenepl->Fill(deltaplanes-0.5,stripPH);
00349   }//loop on strips
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   // loop over strips in current event
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     // 5 plane window
00386     if(ntpStrip->plane >= vtxPlane &&
00387        ntpStrip->plane < vtxPlane + 5){
00388       stpTimes.push_back(ntpStrip->time1);
00389     }
00390   } // end for ns
00391 
00392    if(tfcount) ToFCorrectedEvtTime /= (double)tfcount;
00393    else ToFCorrectedEvtTime=ANtpDefVal::kDouble;
00394 
00395   // calculate strip time median for this event
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   //the event number passed into the method is the current event we care about
00415   evtManip->SetEventInSnarl(event);
00416 
00417   const NtpSREvent* ntpEvent = evtManip->GetEvent();
00418 
00419   // VERTEX TIMING AND DELTA Z CALCULATIONS
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   // Smallest difference in event time to another event
00429   Double_t timeDiff = ANtpDefVal::kFloat;
00430   // Index of the event that had the lowest timeDiff
00431   Int_t closeEventIndex = ANtpDefVal::kInt;
00432   // Smallest different in vtx z to another event
00433   Float_t minDeltaZ = ANtpDefVal::kFloat;
00434   // Minimum time-of-flight-corrected time difference
00435   Double_t mintime = ANtpDefVal::kFloat;
00436  
00437   // Loop over all the other events to find the nearest in time and space
00438   for(int i = 0; i <= ntpManipulator->GetEventArray()->GetLast(); ++i){
00439     // don't want to compare time to itself
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     // time of flight correction for second event
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     // find minimum DeltaZ
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   } // end for i
00473 
00474   eventInfo->minTimeSeparation = timeDiff;
00475   eventInfo->closeTimeEvent = closeEventIndex;
00476   eventInfo->minDeltaZ = minDeltaZ;
00477   eventInfo->medianTime = thisEvtTime;
00478   eventInfo->evtTimeDiff = mintime;
00479   
00480   // ADJACENT STRIP VARIABLES
00481 
00482   // Map from plane number to set of strip numbers in that plane
00483   map<int, set<int> > mplanes;
00484   // Fill map of strips for all other events in this snarl
00485   for(int i = 0; i <= ntpManipulator->GetEventArray()->GetLast(); ++i){
00486     // don't want to compare time to itself
00487     if(i == event) continue;
00488     evtManip->SetEventInSnarl(i);
00489     const NtpSREvent* otherEvent = evtManip->GetEvent();
00490 
00491     // Fill mplanes
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     } // end for ns
00500   } // end for i
00501 
00502   // Number of strips in the event
00503   int nstrips =0;
00504   // Count number of strips with neighbouring strips
00505   int sharedst=0;
00506 
00507   // loop over strips in current event
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     //if no hits on same plane, go to next strip
00519       if(PlanesIter == mplanes.end()) continue;
00520       set<int> stps = PlanesIter->second;
00521       const unsigned int stripn = ntpStrip->strip;
00522       
00523       //when I fill histograms i want only min separation per event
00524       for(set<int>::iterator stpit = stps.begin();
00525           stpit != stps.end();
00526           ++stpit){
00527         // Only neighbouring strips
00528         if(abs(int(stripn)-int(*stpit)) <= 1){
00529           ++sharedst;
00530           break;
00531         }
00532       } // end for stpit
00533     } // end for ns
00534 
00535   if(nstrips) 
00536     eventInfo->sharedStripFraction = float(sharedst)/nstrips;
00537   
00538   // SLICE PH FRACTION
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   //NUMBER OF CONSECUTIVE PLANES
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   //max series of consecutive planes
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   //store the last value it had 
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   //steal contiguous planes variables from NueAna 
00605   //(the implementation is as faithful as possible to the original in NueAna/ShwfitAna.cxx)
00606   //  contPlaneCount                                                                                                   
00607 
00608   //loop over event strips and fill plane histogram ph weighted
00609   //ph is in pes here, in MIP in original nuecode  
00610 
00611   //nophcut
00612   eventInfo->contPlaneCount=contiguousPlanes(ntpEvent,0);
00613   //2.00pe
00614   eventInfo->contPlaneCount050=contiguousPlanes(ntpEvent,2.0);
00615   
00616   // EDGE ACTIVITY VARIABLES
00617 
00618   const double activityTimeStart = thisEvtTime - 4e-8; // +/- 40ns window
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   // loop over strips in snarl
00627   for(int i = 0; i <= fStripArray->GetLast(); ++i){
00628     const NtpSRStrip* strip = (NtpSRStrip*)fStripArray->At(i);
00629 
00630     // consider only strips in time window
00631     if(strip->time1 <= activityTimeStart || strip->time1 >= activityTimeStop)
00632       continue;
00633 
00634     // Only want calorimeter strips
00635     if(strip->plane >= 121) continue;
00636 
00637     if(strip->ph1.pe < 2) continue;
00638 
00639       // look at U planes in calorimeter
00640       if(strip->plane % 2 == 1 && strip->tpos < -0.24){
00641         ++edgeActivityStrips;
00642         edgeActivityPH += strip->ph1.sigcor;
00643       }
00644       // look at opposite edge in U (use 3 strips of fully instr.)
00645       if(strip->plane % 2 == 1 && strip->tpos > 2.27){
00646         ++oppEdgeStrips;
00647         oppEdgePH += strip->ph1.sigcor;
00648       }
00649       // look at V planes in calorimeter
00650       if(strip->plane % 2 == 0 && strip->tpos > 0.24){
00651         ++edgeActivityStrips;
00652         edgeActivityPH += strip->ph1.sigcor;
00653       }
00654       // look at opposite edge in V (use 3 strips of fully instr.)
00655       if(strip->plane % 2 == 0 && strip->tpos <- 2.27){
00656         ++oppEdgeStrips;
00657         oppEdgePH += strip->ph1.sigcor;
00658       }
00659    } // end for i
00660 
00661   eventInfo->edgeActivityPH = edgeActivityPH;
00662   eventInfo->edgeActivityStrips = edgeActivityStrips;
00663   eventInfo->oppEdgePH = oppEdgePH;
00664   eventInfo->oppEdgeStrips = oppEdgeStrips;
00665 
00666   // LARGEST EVENT
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   //make sure to set the event manipulator back to the correct event
00683   ntpManipulator->GetEventManipulator()->SetEventInSnarl(event);
00684 }
00685 
00686 
00687 //----------------------------------------------------------------------
00688 void ANtpInfoObjectFillerNC::InitializekNN(ANtpRecoNtpManipulator *ntpManipulator)
00689 {
00690 
00691   if(!fkNNSet){
00692 
00693     //set the location of the weight files based on the detector
00694     TString baseConf = getenv("SRT_PUBLIC_CONTEXT");
00695     
00696     //check that the public context exists, if not assert false
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     //kreg->Set("PrintConfig", "yes");
00714     kreg->Set("AlgSnarlName", "FillkNN");
00715     kreg->Set("FillkNNKeySignal", int(1));      // for historical reasons muon have type=1                   
00716     //kreg->Set("FillkNNKNeighbor", k);
00717     kreg->Set("FillkNNKNeighbor", 99);
00718     kreg->Set("FillkNNKNeighborMod", int(5));
00719     kreg->Set("FillkNNAddEvent", "no");
00720     //kreg->Set("FillkNNUseTrack", usetrack.c_str());
00721     kreg->Set("FillkNNUseTrack", "no");
00722     kreg->Set("FillkNNTrim", "yes");
00723     kreg->Set("FillkNNTrimDelta", int(1000));
00724     //kreg->Set("FillkNNFilePath", shortknn.c_str());
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     //kreg-> Set("FillkNNKeyTruth", keytrue);
00730     kreg->Set("FillkNNKeyTruth", 19000);
00731     //kreg->Set("FillkNNKeyBase", keybase);
00732     kreg->Set("FillkNNKeyBase", 19400);
00733     //kreg->Set("FillkNNKeyList", keylist.c_str());
00734     kreg->Set("FillkNNKeyList", "19008 19007 19004");
00735     kreg->Set("FillShortEventKeyPass",14001);
00736     kreg->Set("FillShortEventKeyPass2",14001); // no access to 11001
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 //this is the method where you fill all variables related to a track
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   // Estimate dtdz of track, using linear regression of (y = Gx +H)
00791   // <G> = (<xy> - <x><y>) / (<x^2> - <x>^2)
00792   Double_t sumZT(0.); Double_t sumZ2(0.);
00793   Double_t sumT(0.); Double_t sumZ(0.);
00794   Int_t nHits(0); // # of digits at FD, # of strips at ND. whew!
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   //loop over the strips in the record and keep track of how many
00816   //are in each plane
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   }//end loop over event strips to see how many are in each plane
00843 
00844 
00845   // Now loop over all strips in the track
00846   // Count the number of strips in each plane for the track
00847   // Also sum up variables needed to fit t(z)
00848 
00849   for(Int_t i = 0; i < ntpTrack->nstrip; ++i){
00850 
00851     //get the NtpSRStrip object
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   } //end loop over strips in track
00892 
00893   trackInfo->xTalkStrips = uXTalkStripCtr+vXTalkStripCtr;
00894 
00895   //....................................................................
00896   // Calculate dtdz
00897   trackInfo->dtdz = 0;
00898   if (nHits > 1){  // VarZ == 0 if nHits == 1; next conditon would fail
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   //loop over all the planes and if the number of strips in the plane for
00907   //the track is the same as the number for the record, then it is a tracklike plane
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   //does the pathlength of the track make sense?
00939   if(trackInfo->length<50.){
00940 
00941     //loop over all strips in the track
00942     trackInfo->totalStrips = (int)ntpTrack->nstrip;
00943     numDoubleEndedStrips = 0;
00944 
00945     stpCtr = 0;
00946     //ntpStrip = dynamic_cast<NtpSRStrip *>(fStripArray->At(ntpTrack->stp[0]));
00947     for(Int_t ns = 0; ns < trackInfo->totalStrips; ++ns){
00948 
00949       //get the index for this strip
00950       if (ntpTrack->stp[ns] >= 0) index = ntpTrack->stp[ns];
00951       else continue;
00952       //get the NtpSRStrip object
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       } // end if the strip had double sided readout
00967     } // end loop over strips in the track
00968 
00969     trackInfo->twoEndStripFraction = 0.;
00970     if(trackInfo->totalStrips>0)
00971       trackInfo->twoEndStripFraction = (1.*numDoubleEndedStrips)/(1.*trackInfo->totalStrips);
00972 
00973   }//end if reasonable track length
00974 
00975   //fill rustem's pid parameter here. Using new interface.
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 //this is the method where you fill all variables related to a shower
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 //    showerInfo->aShwph_ns        = ntpShower->ph.sigcor;
01046 //    showerInfo->aShwphper_ns        = ntpShower->ph.sigcor/ntpEvent ->ph.sigcor;
01047 //    showerInfo->aShwphperpl_ns      = ntpShower->ph.sigcor/ntpShower->plane.n;
01048 //    showerInfo->aShwphperdig_ns     = ntpShower->ph.sigcor/ntpShower->ndigit;
01049 //    showerInfo->aShwphperstp_ns     = ntpShower->ph.sigcor/ntpShower->nstrip;
01050 
01051 
01052   for(Int_t ns = 0; ns < ntpShower->nstrip; ++ns){
01053 
01054     //get the index for this strip
01055     index = ntpShower->stp[ns];
01056 
01057     //get the NtpSRStrip object
01058     ntpStrip = dynamic_cast<NtpSRStrip *>(fStripArray->At(index));
01059     plane = (int)ntpStrip->plane;
01060     strip = ntpStrip->strip;
01061 
01062     //......................................................................
01063     // -=NEW VARIABLES=-
01064     // use all strips for now
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     }//end if far detector
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   }//end loop over strips
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   // -=NEW VARIABLES=-
01122   // fills the variables emFrac and phKurtosis
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     // Kurtosis of ph distribution as defined in Numerical Recipes
01139     showerInfo->phKurtosis = fourthMoment/(nstrip*sigma*sigma)-3.0;
01140 
01141     Int_t emLikeStrips(0);
01142     Bool_t cluFailure(false); // check for very rare cluster failure
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 //this is the method where you fill all variables related to a showers and tracks
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       //get the index for this strip
01191       indexShowerStrip = ntpShower->stp[ns];
01192 
01193       //loop over the track strips
01194       for(Int_t nt = 0; nt < ntpTrack->nstrip; ++nt){
01195 
01196         indexTrackStrip = ntpTrack->stp[nt];
01197 
01198         //check to see if the shower index matches the track index, if so,
01199         //increment the counter for the number of track strips in the shower
01200         if(indexTrackStrip == indexShowerStrip) ++eventInfo->trackStripsInShower;
01201 
01202       }//end loop over track strips
01203 
01204     }//end loop over shower strips
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   }//end if shower and track are both present and track length < 50. m
01212 
01213   //fill niki's variables
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   // loop over strips to compute ph3 ph6 phlast phcommon
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     // shower strips
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     // tracks strips
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    } // end of if stp_index != -1
01311 
01312   } // end of looping over event strips....
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 //this is the method where you fill all variables related to an event
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   //add up the pulse height in each plane to find the maximum in the event
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   // -=NEW VARIABLES=-
01359 
01360   // PLANE ADC (`GLOBAL') VARIABLES
01361   // First make summary of plane ADC
01362   const UInt_t nPlanes(eventInfo->planes);
01363   vector<Float_t> planeAdc(nPlanes+12,0);
01364   // pad vector with 12 zeros. these are needed later
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   // Also calculate some simple event duration type variables
01369   Double_t t_1stStp(10.); //time of first strip in event
01370   Double_t t_lastStp(0.); //  & of last
01371   Double_t t_sumStp(0.);  // sum of times
01372   Double_t t_sum2Stp(0.); // sum of squares of times
01373   Int_t nStpHit( ntpEvent->nstrip ); // number of double ended strips
01374 //......................................................................
01375 
01376   for(Int_t ns = 0; ns < ntpEvent->nstrip; ++ns){
01377 
01378     //get the current strip
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     // -=NEW VARIABLES=-
01402 
01403     // Plane ADC vars
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     // event duration variables
01411     Double_t t_stpHitA = TMath::Max(ntpStrip->time0 , ntpStrip->time1);
01412     Double_t t_stpHitB = TMath::Min(ntpStrip->time0 , ntpStrip->time1);
01413     //Earlier hit may be -ive, meaning a single ended hit
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   }//end loop over event strips
01429 
01430   //......................................................................
01431   // -=NEW VARIABLES=-
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   //find the rms of the ph per plane for the event and the
01447   //plane with the max ph.  also find the maximum signal in windows of
01448   //3,6,9,12 planes from the first plane
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     } //end loop over planes in event
01478 
01479   }//end loop over window size
01480 
01481   eventInfo->pulseHeightRms = planePHDist.GetRMS();
01482 
01483   //......................................................................
01484   // -=NEW VARIABLES=-
01485 
01486   //Calculate pulse height in plane groups 1 to 3 and 4 to 6
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   // Image recognition discriminants (KEA)
01498   //FillVHSEvtInfo(ntpManipulator, eventInfo);
01499 
01500   return;
01501 }
01502 
01503 //----------------------------------------------------------------------
01504 //this is the method where you fill all variables related to the MC truth
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   //fill the plane to pmt map first.  the idea is that there are at most 2 PMT's per plane
01527   //in either detector.  put one NtpSRStrip->pmtindex0(1) into the fPlaneToPMTMapEast(West)
01528   //array.  near detector only has pmt's on the west side.  call that pmt 0 for the plane.
01529   //then when filling the signal arrays if the pmtindex0(1) for a strip is not the one in the
01530   //array, it is pmt 1 for the plane.
01531   for(Int_t ns = 0; ns < ntpEvent->nstrip; ++ns){
01532 
01533     //get the NtpSRStrip object
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   }//end loop over strips to fill plane to pmt maps
01550 
01551   MSG("ANtpInfoObjectFillerNC", Msg::kDebug)
01552     << "finish filling plane to pmt maps" << endl;
01553 
01554   //loop over all the strips and put the signal size in the arrays
01555   for(Int_t ns = 0; ns < ntpEvent->nstrip; ++ns){
01556 
01557     //get the NtpSRStrip object
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     //got the strip, now fill the plane pixel info
01569     if(fDetector == Detector::kFar){
01570 
01571       //east first
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       //now west
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     }//end if far detector;
01587     else if(fDetector == Detector::kNear){
01588 
01589       //figure out what type of plane you are
01590       //planeview = 2 --> U, 3 --> V
01591       //map types 2 and 3 are for full planes, 1 is for partial
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       }//end if full
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       }//end if partial
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     }//end if near detector
01636   }//end loop over strips
01637 
01638   return;
01639 }
01640 
01641 //----------------------------------------------------------------------
01642 void ANtpInfoObjectFillerNC::FindEarlyActivityWeight(NtpSREvent *ntpEvent,
01643                                                      ANtpEventInfoNC *eventInfo)
01644 {
01645 
01646   //find the time of the earliest strip in the event
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 //   for(int i = 0; i < ntpEvent->nstrip; ++i){
01662 //     strip = dynamic_cast<NtpSRStrip *>(fStripArray->At(ntpEvent->stp[i]));
01663 
01664 //     MSG("ANtpInfoObjectFillerNC", Msg::kInfo)
01665 //       << fHeaderInfo->snarl << " " << fHeaderInfo->events
01666 //       << " " << eventInfo->event << " "
01667 //       <<strip->plane << " " << earliestEventTime
01668 //       << " " << 1.e9*(strip->time1-earliestEventTime) << endl;
01669 //   }
01670 
01671   //now find the weigthed sum of ADC for the early strips - make sure the early
01672   //strips come from the same pmts as the event strips.
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     //work in ns
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   }//end if far detector
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     //get the NtpSRStrip object
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   }//end loop over strips
01824   return;
01825 }
01826