loon: NtpEMModule.cxx Source File

00001 
00002 //
00003 // NtpEMModule.cxx
00004 //
00005 // A JobControl Module for filling an NtpEMRecord
00006 //
00008 #include <iostream>
00009 using namespace std;
00010 
00011 #include "TClonesArray.h"
00012 #include "CandNtupleEM/Module/NtpEMModule.h"
00013 #include "CandNtupleEM/NtpEMRecord.h"
00014 #include "CandNtupleEM/NtpEMShower.h"
00015 #include "CandNtupleEM/NtpEMVertex.h"
00016 #include "CandNtupleEM/NtpEMSummary.h"
00017 #include "CandNtupleSR/NtpSRStrip.h"
00018 #include "CandNtupleSR/NtpSRSlice.h"
00019 #include "CandNtupleSR/NtpSRDate.h"
00020 #include "CandNtupleSR/NtpSRPlane.h"
00021 #include "MessageService/MsgService.h"
00022 #include "JobControl/JobCModuleRegistry.h"
00023 #include "JobControl/JobCommand.h"
00024 #include "MinosObjectMap/MomNavigator.h"
00025 #include "RawData/RawDigit.h"
00026 #include "RawData/RawRecord.h"
00027 #include "RawData/RawDaqSnarlHeader.h"
00028 #include "RawData/RawDigitDataBlock.h"
00029 #include "Record/RecCandHeader.h"
00030 #include "CandData/CandRecord.h"
00031 #include "CandData/CandHeader.h"
00032 #include "RecoBase/CandStripListHandle.h"
00033 #include "RecoBase/CandStripHandle.h"
00034 #include "RecoBase/CandSliceListHandle.h"
00035 #include "RecoBase/CandSliceHandle.h"
00036 #include "RecoBase/CandShowerListHandle.h"
00037 #include "RecoBase/CandShowerHandle.h"
00038 #include "RecoBase/PropagationVelocity.h"
00039 #include "CandShowerEM/CandShowerEMHandle.h"
00040 #include "CandShowerEM/CandShowerEMListHandle.h"
00041 #include "CandFitShowerEM/CandFitShowerEMHandle.h"
00042 #include "CandDigit/CandDeMuxDigitListHandle.h"
00043 #include "UgliGeometry/UgliGeomHandle.h"
00044 #include "Plex/PlexHandle.h"
00045 #include "Plex/PlexVetoShieldHack.h"
00046 
00047 #include <cassert>
00048 
00049 ClassImp(NtpEMModule)
00050 
00051 //   Definition of static data members
00052 //   *********************************
00053 
00054 CVSID("$Id: NtpEMModule.cxx,v 1.6 2005/05/18 03:50:53 schubert Exp $");
00055 JOBMODULE(NtpEMModule, "NtpEMModule",
00056          "A module for filling EM ntuple records.");
00057 
00058 const Double_t NtpEMModule::kCos45 = 0.70710678;  //used to convert u,v to x,y
00059 
00060 //......................................................................
00061 
00062 const Registry& NtpEMModule::DefaultConfig() const {
00063 
00064   MSG("NtpEM",Msg::kDebug) << 
00065     "NtpEMModule::DefaultConfig" << endl;
00066 
00067   static Registry r; 
00068   std::string name = this->JobCModule::GetName();
00069   name += ".config.default";
00070   r.SetName(name.c_str());
00071 
00072   r.UnLockValues();
00073   r.Set("WriteStrip",false);
00074   r.Set("WriteSlice",false);
00075   r.LockValues();
00076   
00077   return r;
00078 }
00079 
00080 //......................................................................
00081 
00082 void NtpEMModule::Config(const Registry& r) {
00083   
00084   MSG("NtpEM",Msg::kDebug) << "NtpEMModule::Config" << endl;
00085   
00086   Int_t tmpi;
00087 
00088   if ( r.Get("WriteStrip",tmpi) ) {
00089     fWriteStrip = tmpi;
00090   }
00091   if ( r.Get("WriteSlice",tmpi) ) {
00092     fWriteSlice = tmpi;
00093   }
00094   
00095 }
00096 
00097 //......................................................................
00098 
00099 JobCResult NtpEMModule::Reco(MomNavigator *mom) {
00100   
00101   JobCResult result(JobCResult::kPassed);  
00102   MSG("NtpEM",Msg::kDebug) << "NtpEMModule::Reco" << endl;
00103 
00104   // Reset maps used to associate uid of reconstructed object with array index
00105   fStripUidMap.clear();
00106   fSliceUidMap.clear();
00107   fShowerUidMap.clear();
00108   
00109   // Check that mom exists.
00110   assert(mom);
00111 
00112   const RawRecord* rawrec = 
00113    dynamic_cast<const RawRecord*>(mom->GetFragment("RawRecord","","DaqSnarl"));
00114   if (!rawrec) {
00115     MSG("NtpEM",Msg::kWarning) << "No DaqSnarl RawRecord in Mom" << endl;
00116     result.SetWarning().SetFailed();
00117     return result;
00118   }
00119 
00120   const CandRecord* cndrec = dynamic_cast<const CandRecord*>
00121           (mom->GetFragment("CandRecord","PrimaryCandidateRecord")); 
00122   if (!cndrec) {
00123     MSG("NtpEM",Msg::kWarning) << "No PrimaryCandidateRecord in Mom" << endl;
00124     result.SetWarning().SetFailed();
00125     return result;
00126   }
00127 
00128   // Extract header from CandRecord and use this to create RecCandHeader
00129   // and NtpEMRecord.
00130   const CandHeader* cndhdr = cndrec -> GetCandHeader();
00131   const RawDaqSnarlHeader* rawhdr = dynamic_cast<const RawDaqSnarlHeader*>
00132                                     (rawrec -> GetRawHeader());
00133   RecCandHeader ntphdr(rawhdr->GetVldContext(),rawhdr->GetRun(),
00134          rawhdr->GetSubRun(),rawhdr->GetRunType(),rawhdr->GetErrorCode(),
00135          rawhdr->GetSnarl(),rawhdr->GetTrigSrc(),rawhdr->GetTimeFrameNum(),
00136          rawhdr->GetRemoteSpillType(),cndhdr->GetEvent());
00137   NtpEMRecord* ntprec = new NtpEMRecord(ntphdr);  
00138 
00139   this -> FillNtpStrip(ntprec,cndrec);
00140   this -> FillNtpSlice(ntprec,cndrec);
00141   this -> FillNtpShower(ntprec,cndrec);  
00142   this -> FillNtpSummary(ntprec,cndrec);
00143 
00144   mom -> AdoptFragment(ntprec);  // pass record to mom to own
00145 
00146   return result;
00147 
00148 }
00149 
00150 void NtpEMModule::FillNtpStrip(NtpEMRecord* ntprec,const CandRecord* cndrec) {
00151   
00152   MSG("NtpEM",Msg::kVerbose) << "NtpEMModule::FillNtpStrip" << endl;
00153 
00154   const CandStripListHandle *striplisthandle 
00155    = dynamic_cast <const CandStripListHandle*> 
00156      (cndrec->FindCandHandle("CandStripListHandle"));
00157   if ( !striplisthandle ) return; // no strips => done
00158   
00159   Int_t nstrip = 0;
00160   TIter stripItr(striplisthandle->GetDaughterIterator());
00161   TClonesArray& striparray = *(ntprec->stp);
00162   while ( CandStripHandle* strip=dynamic_cast<CandStripHandle*>(stripItr())) {
00163     Int_t uid = strip->GetUidInt();
00164     fStripUidMap.insert(std::make_pair(uid,nstrip));    
00165     if(fWriteStrip){
00166       // Uses new with placement to efficiently create strip ntp
00167       NtpSRStrip* ntpstrip = new((striparray)[nstrip++])NtpSRStrip();
00168       // Transport information from CandStrip to strip ntp
00169       ntpstrip->index = nstrip-1;
00170       ntpstrip->planeview = strip->GetPlaneView(); // plane view
00171       ntpstrip->ndigit = strip->GetNDigit();
00172       ntpstrip->demuxveto = strip->GetDemuxVetoFlag();
00173       ntpstrip->strip = strip->GetStrip(); // strip number
00174       ntpstrip->plane = strip->GetPlane(); // plane number
00175       ntpstrip->tpos = strip->GetTPos();  // tpos
00176       ntpstrip->z = strip->GetZPos();  // zpos
00177       
00178       // Raw channel id of first digit associated with each end
00179       bool negEndDone = false;
00180       bool posEndDone = false;
00181       TIter digitItr(strip -> GetDaughterIterator());
00182       while (CandDigitHandle* digit=dynamic_cast<CandDigitHandle*>(digitItr())) {
00183         const RawChannelId& rawch = digit->GetChannelId();
00184         Int_t pmtindex = rawch.GetCrate()*108 + rawch.GetVarcId()*36
00185           +rawch.GetVmm()*6 + rawch.GetVaAdcSel()*3+rawch.GetVaChip();
00186         if ( !negEndDone && 
00187              digit -> GetPlexSEIdAltL().GetEnd()==StripEnd::kNegative) { 
00188           ntpstrip->SetPmtIndex(pmtindex,0); negEndDone = true;  
00189         }
00190         else if( !posEndDone && 
00191                  digit -> GetPlexSEIdAltL().GetEnd()==StripEnd::kPositive) {
00192           ntpstrip->SetPmtIndex(pmtindex,1); posEndDone = true;
00193         }  
00194       }
00195       
00196       // Strip end dependent quantities
00197       for (UInt_t i = 0; i < 2; i++ ) {
00198         StripEnd::EStripEnd stripend = StripEnd::kNegative;
00199         if (i == 1) stripend = StripEnd::kPositive;
00200         if ( strip->GetNDigit(stripend) > 0 ) {
00201           NtpSRPulseHeight ph;
00202           ph.raw = strip->GetCharge(CalDigitType::kNone,stripend);
00203           ph.siglin = strip->GetCharge(CalDigitType::kSigLin,stripend);
00204           ph.sigcor = strip->GetCharge(CalDigitType::kSigCorr,stripend);
00205           ph.pe = strip->GetCharge(CalDigitType::kPE,stripend);
00206           ntpstrip->SetPh(ph,i); 
00207           ntpstrip->SetTime(strip->GetTime(stripend),i);
00208         }
00209       }
00210     }
00211     else nstrip++;
00212   } // done with all strips
00213 
00214   return;
00215 
00216 }
00217 
00218 void NtpEMModule::FillNtpSlice(NtpEMRecord* ntprec,const CandRecord* cndrec) {
00219 
00220   MSG("NtpEM",Msg::kVerbose) << "NtpEMModule::FillNtpSlice" << endl;
00221   
00222   const CandSliceListHandle *slicelisthandle 
00223     = dynamic_cast <const CandSliceListHandle*> 
00224     (cndrec -> FindCandHandle("CandSliceListHandle"));
00225   if ( !slicelisthandle ) return; // all done
00226   
00227   Int_t nslice = 0;
00228   TIter sliceItr(slicelisthandle->GetDaughterIterator());
00229   TClonesArray& slicearray = *(ntprec->slc);
00230   while ( CandSliceHandle* slice=dynamic_cast<CandSliceHandle*>(sliceItr())) {
00231     Int_t uid = slice->GetUidInt();
00232     fSliceUidMap.insert(std::make_pair(uid,nslice));
00233     if(fWriteSlice){
00234       // Uses new with placement to efficiently create slice ntp    
00235       NtpSRSlice* ntpslice = new(slicearray[nslice++])
00236         NtpSRSlice(slice->GetNStrip());
00237       ntpslice->index  = nslice-1; // index is number of slices - 1
00238       ntpslice->ndigit = slice->GetNDigit();
00239       // Fill indices of associated strips in slice tree
00240       TIter slicestripItr(slice->GetDaughterIterator());
00241       Int_t nslicestrip = 0;
00242       while ( CandStripHandle *slicestrip = dynamic_cast<CandStripHandle*>
00243               (slicestripItr())) {
00244         Int_t uid = slicestrip->GetUidInt();
00245         Int_t stripindex = fStripUidMap[uid];
00246         ntpslice->AddStripAt(stripindex,nslicestrip); // add index to strip
00247         nslicestrip++;
00248       }
00249       
00250       // Set summed charge in slice
00251       ntpslice->ph.raw = slice->GetCharge(CalDigitType::kNone);
00252       ntpslice->ph.siglin = slice->GetCharge(CalDigitType::kSigLin);
00253       ntpslice->ph.sigcor = slice->GetCharge(CalDigitType::kSigCorr);
00254       ntpslice->ph.pe = slice->GetCharge(CalDigitType::kPE);
00255       
00256       // Set range of planes included in slice
00257       ntpslice->plane.n   =slice->GetNPlane();
00258       ntpslice->plane.beg = slice->GetBegPlane();
00259       ntpslice->plane.end = slice->GetEndPlane();
00260       ntpslice->plane.nu   = slice->GetNPlane(PlaneView::kU);
00261       ntpslice->plane.begu = slice->GetBegPlane(PlaneView::kU);
00262       ntpslice->plane.endu = slice->GetEndPlane(PlaneView::kU);
00263       ntpslice->plane.nv   = slice->GetNPlane(PlaneView::kV);
00264       ntpslice->plane.begv = slice->GetBegPlane(PlaneView::kV);
00265       ntpslice->plane.endv = slice->GetEndPlane(PlaneView::kV);
00266       
00267     }
00268   }
00269   return;
00270 }
00271 
00272 void NtpEMModule::FillNtpShower(NtpEMRecord* ntprec,const CandRecord* cndrec) {
00273 
00274   MSG("NtpEM",Msg::kVerbose) << "NtpEMModule::FillNtpShower" << endl;
00275 
00276   const CandShowerEMListHandle *showerlisthandle 
00277    = dynamic_cast <const CandShowerEMListHandle*> 
00278      (cndrec -> FindCandHandle("CandShowerEMListHandle"));
00279   if ( !showerlisthandle ) return; // all done
00280 
00281   Int_t nshower = 0;
00282   TIter showerItr(showerlisthandle->GetDaughterIterator());
00283   TClonesArray& showerarray = *(ntprec->emshw);
00284   while (CandShowerEMHandle* shower = 
00285          dynamic_cast<CandShowerEMHandle*> (showerItr())) {
00286     
00287     Int_t uid = shower->GetUidInt();
00288     fShowerUidMap.insert(std::make_pair(uid,nshower));
00289     
00290     // Uses new with placement to efficiently create slice ntp    
00291     NtpEMShower* ntpshower = new(showerarray[nshower++])
00292       NtpEMShower(shower->GetNStrip());
00293     // Fill indices of associated strips in shower tree
00294     ntpshower->index = nshower - 1;
00295     ntpshower->ndigit = shower->GetNDigit();
00296     
00297     TIter showerstripItr(shower->GetDaughterIterator());
00298     Int_t nshowerstrip = 0;
00299     
00300     while (CandStripHandle *showerstrip = 
00301            dynamic_cast<CandStripHandle*> (showerstripItr())) {
00302  
00303       Int_t uid = showerstrip->GetUidInt();
00304       Int_t stripindex = fStripUidMap[uid];
00305       
00306       // add index to strip
00307       ntpshower->AddStripAt(stripindex,nshowerstrip); 
00308       nshowerstrip++;
00309       
00310     }
00311 
00312     // Set range of planes included in slice
00313     ntpshower->plane.n = shower->GetNPlane();
00314     ntpshower->plane.nu = shower->GetNPlane(PlaneView::kU);
00315     ntpshower->plane.nv = shower->GetNPlane(PlaneView::kV);
00316     ntpshower->plane.beg = shower->GetBegPlane();
00317     ntpshower->plane.begu = shower->GetBegPlane(PlaneView::kU);
00318     ntpshower->plane.begv = shower->GetBegPlane(PlaneView::kV);
00319     ntpshower->plane.end = shower->GetEndPlane();
00320     ntpshower->plane.endu = shower->GetEndPlane(PlaneView::kU);
00321     ntpshower->plane.endv = shower->GetEndPlane(PlaneView::kV);
00322  
00323     // Set summed charge in shower
00324     ntpshower->ph.raw = shower->GetCharge(CalStripType::kNone);
00325     ntpshower->ph.siglin = shower->GetCharge(CalStripType::kSigLin);
00326     ntpshower->ph.sigcor = shower->GetCharge(CalStripType::kSigCorr);
00327     ntpshower->ph.pe = shower->GetCharge(CalStripType::kPE);
00328     ntpshower->ph.sigmap = shower->GetCharge(CalStripType::kSigMapped);
00329     ntpshower->ph.mip = shower->GetCharge(CalStripType::kMIP);
00330     ntpshower->ph.gev = shower->GetCharge(CalStripType::kGeV);
00331  
00332     // Set vertex of shower
00333     NtpEMVertex& vtx = ntpshower->vtx;
00334     vtx.u = shower->GetVtxU();
00335     vtx.v = shower->GetVtxV();
00336     vtx.x = kCos45*(vtx.u-vtx.v);
00337     vtx.y = kCos45*(vtx.u+vtx.v);
00338     vtx.z = shower->GetVtxZ();
00339     vtx.t = shower->GetVtxT();
00340     vtx.plane = shower->GetVtxPlane();
00341     vtx.dcosu = shower->GetDirCosU();
00342     vtx.dcosv = shower->GetDirCosV();
00343     vtx.dcosx = kCos45*(vtx.dcosu-vtx.dcosv);
00344     vtx.dcosy = kCos45*(vtx.dcosu+vtx.dcosv);
00345     vtx.dcosz = shower->GetDirCosZ(); 
00346 
00347     //set other things specific to EMShower:
00348     const Double_t *temp = shower->GetEigenVectors();
00349     for(int i=0;i<8;i++) ntpshower->evector[i] = temp[i];
00350     temp = shower->GetEigenValues();
00351     for(int i=0;i<4;i++) ntpshower->evalue[i] = temp[i];
00352     temp = shower->GetAvgDev();
00353     for(int i=0;i<4;i++) ntpshower->avgdev[i] = temp[i];
00354     temp = shower->GetOutPH();
00355     for(int i=0;i<5;i++) ntpshower->outph[i] = temp[i];
00356     ntpshower->status = shower->GetShwStatus();
00357    
00358     if(CandFitShowerEMHandle* fitshower = 
00359        dynamic_cast<CandFitShowerEMHandle*> (shower)) {
00360       ntpshower->fitinfo.pass = fitshower->GetPass();
00361       ntpshower->fitinfo.momentum = fitshower->GetMomentum();
00362       ntpshower->fitinfo.chi2 = fitshower->GetChi2();
00363       ntpshower->fitinfo.ndof = fitshower->GetNDOF();
00364       ntpshower->fitinfo.cputime = fitshower->GetCPUTime();
00365       ntpshower->fitinfo.niterate = fitshower->GetNIterate();
00366     }
00367   }
00368 
00369   return;
00370 }
00371 
00372 void NtpEMModule::FillNtpSummary(NtpEMRecord* ntprec,
00373                                  const CandRecord* cndrec) {
00374 
00375   MSG("NtpEM",Msg::kVerbose) << "NtpEMModule::FillNtpSummary" << endl;
00376   const CandDigitListHandle *dlh = dynamic_cast <const CandDigitListHandle*> 
00377      (cndrec->FindCandHandle("CandDigitListHandle"));
00378   if ( !dlh ) return; // no digits => done
00379 
00380   NtpEMSummary& ntpsummary = ntprec->emhdr;
00381   ntpsummary.nshower = fShowerUidMap.size();
00382   ntpsummary.nstrip  = fStripUidMap.size();
00383   ntpsummary.nslice  = fSliceUidMap.size();
00384   
00385   ntpsummary.trigtime = dlh->GetAbsTime();
00386   // Fill the NtpEMDate portion of the NtpEMSummary object
00387   const VldContext& vld = *(cndrec->GetVldContext());
00388   UInt_t year,month,day,hour,minute,sec;
00389   vld.GetTimeStamp().GetDate(kTRUE,0,&year,&month,&day);
00390   vld.GetTimeStamp().GetTime(kTRUE,0,&hour,&minute,&sec);
00391   ntpsummary.date.year   = (UShort_t)year;  
00392   ntpsummary.date.month  = (Char_t)month;   
00393   ntpsummary.date.day    = (Char_t)day;     
00394   ntpsummary.date.hour   = (Char_t)hour;  
00395   ntpsummary.date.minute = (Char_t)minute; 
00396   ntpsummary.date.sec    =  (Double_t)sec; 
00397   ntpsummary.date.sec   += ntpsummary.trigtime; 
00398   ntpsummary.date.utc = vld.GetTimeStamp().GetSec();
00399   
00400   Int_t minplaneall  = -1;
00401   Int_t minplaneallu = -1;
00402   Int_t minplaneallv = -1;
00403   Int_t maxplaneall  = -1;
00404   Int_t maxplaneallu = -1;
00405   Int_t maxplaneallv = -1;
00406   std::map<Int_t,Bool_t> planeoccupancyall; 
00407   std::map<Int_t,Bool_t> planeoccupancyallu;
00408   std::map<Int_t,Bool_t> planeoccupancyallv;
00409   Float_t planepe[1000][2] = {{0},{0}}; // initializes all members to 0
00410   
00411   // Loop over all digits
00412   TIter digitItr(dlh -> GetDaughterIterator());
00413   while (CandDigitHandle* digit = dynamic_cast<CandDigitHandle*>(digitItr())) {
00414     // Calculate the summed pulse height of all digits (non-shield) in 
00415     // the entire event
00416     PlexSEIdAltL pseid(digit->GetPlexSEIdAltL());
00417     if (pseid.IsVetoShield()) continue;
00418     ntpsummary.ndigit++;
00419     
00420     ntpsummary.ph.raw += digit->GetCharge(CalDigitType::kNone);
00421     Float_t calcharge[3] = {0};
00422     if ( pseid.GetDemuxVetoFlag() == 0 ) { 
00423       // demux successful
00424       calcharge[0] = digit->GetCharge(CalDigitType::kSigLin);
00425       calcharge[1] = digit->GetCharge(CalDigitType::kSigCorr);
00426       calcharge[2] = digit->GetCharge(CalDigitType::kPE);
00427     }
00428     else if ( pseid.GetSize() > 0 ) {
00429       // if it wasn't demuxed, then simply use first entry
00430       calcharge[0] = pseid[0].GetSigLin();
00431       calcharge[1] = pseid[0].GetSigCorr();
00432       calcharge[2] = pseid[0].GetPE();
00433     }
00434     ntpsummary.ph.siglin += calcharge[0];
00435     ntpsummary.ph.sigcor += calcharge[1];
00436     ntpsummary.ph.pe     += calcharge[2];
00437     
00438     // Now determine the range of planes.  In the first case, the range of
00439     // planes is determined as the min/max plane with any digit.  
00440     Int_t iplane = pseid.GetPlane();
00441     if ( iplane < 0 || iplane >= 1000 ) continue;
00442     
00443     if ( minplaneall < 0 || iplane < minplaneall ) minplaneall = iplane;
00444     if ( maxplaneall < 0 || iplane > maxplaneall ) maxplaneall = iplane;
00445     planeoccupancyall[iplane] = kTRUE; // at least one digit on this plane
00446     
00447     switch (pseid.GetPlaneView()) {
00448     case PlaneView::kU:
00449       if (minplaneallu < 0 || iplane < minplaneallu) minplaneallu = iplane;
00450       if (maxplaneallu < 0 || iplane > maxplaneallu) maxplaneallu = iplane;
00451       planeoccupancyallu[iplane] = kTRUE; // at least one digit on this u-plane
00452       break;
00453     case PlaneView::kV:
00454       if (minplaneallv < 0 || iplane < minplaneallv) minplaneallv = iplane;
00455       if (maxplaneallv < 0 || iplane > maxplaneallv) maxplaneallv = iplane;
00456       planeoccupancyallv[iplane] = kTRUE; // at least one digit on this v-plane
00457       break;
00458     default:
00459       break;
00460     }
00461     
00462     // In the second case, the range of plane requires determining which planes
00463     // have a summed ph (over both readout ends) of > 3 pe. Store the plane
00464     // pe sum now and use it below.
00465     switch (pseid.GetEnd()) {
00466     case StripEnd::kNegative:
00467       planepe[iplane][0] += calcharge[2];  // CalDigitType::kPE
00468       break;
00469     case StripEnd::kPositive:
00470       planepe[iplane][1] += calcharge[2];  // CalDigitType::kPE
00471       break;
00472     default:
00473       break;
00474     }
00475   } // end of digit while loop
00476   
00477   ntpsummary.planeall.beg  = minplaneall;
00478   ntpsummary.planeall.end  = maxplaneall;
00479   ntpsummary.planeall.begu = minplaneallu;
00480   ntpsummary.planeall.endu = maxplaneallu;
00481   ntpsummary.planeall.begv = minplaneallv;
00482   ntpsummary.planeall.endv = maxplaneallv;
00483   std::map<Int_t,Bool_t>::iterator iter;
00484   for ( Int_t iplane = minplaneall; iplane <= maxplaneall; iplane ++ ) {
00485     iter = planeoccupancyall.find(iplane); // at least one digit
00486     if ( iter != planeoccupancyall.end() ) ntpsummary.planeall.n++;
00487     iter = planeoccupancyallu.find(iplane); // at least one digit u
00488     if ( iter != planeoccupancyallu.end() ) ntpsummary.planeall.nu++;
00489     iter = planeoccupancyallv.find(iplane); // at least one digit v
00490     if ( iter != planeoccupancyallv.end() ) ntpsummary.planeall.nv++;
00491   }
00492   
00493   // plane.beg/end is determined by first testing for 4 contiguous planes
00494   // with a summed ph > 3 pe across both ends. The minimum plane of the first 
00495   // such group and the maximum plane of the last such group form plane.beg 
00496   // and plane.end respectively.
00497   Int_t minplane  = -1;
00498   Int_t maxplane  = -1;
00499   bool found(0);
00500   for (Int_t iplane=minplaneall;iplane<=maxplaneall-3 && !found; iplane++) {
00501     if ( !found && (planepe[iplane][0] + planepe[iplane][1]) > 3.  
00502          && (planepe[iplane+1][0] + planepe[iplane+1][1]) > 3. 
00503          && (planepe[iplane+2][0] + planepe[iplane+2][1]) > 3. 
00504          && (planepe[iplane+3][0] + planepe[iplane+3][1]) > 3. ) {
00505       found = 1;
00506       minplane = iplane;
00507     }
00508   }
00509   found = 0;
00510   for (Int_t iplane=maxplaneall;iplane>=minplaneall+3 && !found; iplane--) {
00511     if ( !found && (planepe[iplane][0] + planepe[iplane][1]) > 3.  
00512          && (planepe[iplane-1][0] + planepe[iplane-1][1]) > 3. 
00513          && (planepe[iplane-2][0] + planepe[iplane-2][1]) > 3. 
00514          && (planepe[iplane-3][0] + planepe[iplane-3][1]) > 3. ) {
00515       found = 1;
00516       maxplane = iplane;
00517     }
00518   }
00519     
00520   Int_t minplaneu = -1;
00521   Int_t minplanev = -1;
00522   Int_t maxplaneu = -1;
00523   Int_t maxplanev = -1;
00524   for (Int_t iplane = minplane; iplane <= maxplane; iplane++ ) {
00525     if ( iplane >= 0 ) {
00526       if ( planepe[iplane][0] + planepe[iplane][1] > 0. ) {
00527         ntpsummary.plane.n++;  // non-zero readout at either end
00528         iter = planeoccupancyallu.find(iplane);
00529         if ( iter != planeoccupancyallu.end() ) {
00530           ntpsummary.plane.nu++; // non-zero readout on u plane
00531           if ( minplaneu < 0 || iplane < minplaneu ) minplaneu = iplane;
00532           if ( maxplaneu < 0 || iplane > maxplaneu ) maxplaneu = iplane;
00533         }
00534         iter = planeoccupancyallv.find(iplane);
00535         if ( iter != planeoccupancyallv.end() ) {
00536           ntpsummary.plane.nv++;
00537           if ( minplanev < 0 || iplane < minplanev ) minplanev = iplane;
00538           if ( maxplanev < 0 || iplane > maxplanev ) maxplanev = iplane;
00539         }
00540       }
00541     }
00542   }
00543   
00544   ntpsummary.plane.beg = minplane;
00545   ntpsummary.plane.end = maxplane;
00546   ntpsummary.plane.begu = minplaneu;
00547   ntpsummary.plane.endu = maxplaneu;
00548   ntpsummary.plane.begv = minplanev;
00549   ntpsummary.plane.endv = maxplanev;
00550   
00551   return;
00552 
00553 }
00554 
00555 
00556 
00557 
00558 
00559