VHS Namespace Reference

Enumerations

enum  evtType {
  vhsNC = 0, vhsCCe = 1, vhsCCmu = 2, vhsCCtau = 3,
  vhsUnknown = 4
}

Functions

bool DrawEvent (int nEntry, const char *treeName, int nPlanes, int nStrips)
void FillDiscriminants (NtpSREvent *evt, TClonesArray *stp, std::vector< double > avgNC, std::vector< double > avgCCe, std::vector< double > avgCCmu, std::vector< double > avgCCtau, std::vector< double > varNC, std::vector< double > varCCe, std::vector< double > varCCmu, std::vector< double > varCCtau, std::vector< double > pNC, std::vector< double > pCCe, std::vector< double > pCCmu, std::vector< double > pCCtau, const int nPlanes, const int nStrips, const bool bUnit, VHSevent *&vhsevt)
std::vector< double > FindMedian (std::vector< std::vector< double > > pts, std::vector< double > avg, bool bUnit=true, bool bVerbose=true)
std::vector< double > FullVec (std::vector< double > image, std::vector< int > index, int nPlanes, int nStrips)
double GetDistance (std::vector< double > vec0, std::vector< double > vec1)
VHS::evtType GetEvtType (int inu, int iaction)
void GetImage (int *index, int nstp, TClonesArray *stp, int nPlanes, int nStrips, std::vector< double > &image, std::vector< int > &vecInd, std::vector< double > &theta)
double GetLL (std::vector< double > fullImage, std::vector< double > pHit, std::vector< double > avg, std::vector< double > var)
int GetPlane (int vecInd, int nPlanes)
int GetStrip (int vecInd, int nPlanes)
void GetThetaAxis (TClonesArray *stp, std::vector< int > index, double &theta, std::vector< int > &center)
std::vector< std::string > ParseNmList (const char *cstr)
void ReadFile (TFile *inFile, std::vector< double > &avgNC, std::vector< double > &avgCCe, std::vector< double > &avgCCmu, std::vector< double > &avgCCtau, std::vector< double > &varNC, std::vector< double > &varCCe, std::vector< double > &varCCmu, std::vector< double > &varCCtau, std::vector< double > &pNChit, std::vector< double > &pCCehit, std::vector< double > &pCCmuhit, std::vector< double > &pCCtauhit, const int nPlanes)
void RotatePixel (int &plane, int &strip, const int avgPlane, const int avgStrip, const double theta)
void SeparateViews (int *index, int nstp, TClonesArray *stp, std::vector< int > &ustp, std::vector< int > &vstp)
std::vector< VHSevent * > Skim (NtpStRecord *ntpst, std::vector< double > avgNC, std::vector< double > avgCCe, std::vector< double > avgCCmu, std::vector< double > avgCCtau, std::vector< double > varNC, std::vector< double > varCCe, std::vector< double > varCCmu, std::vector< double > varCCtau, std::vector< double > pNC, std::vector< double > pCCe, std::vector< double > pCCmu, std::vector< double > pCCtau, const int nPlanes=20, const int nStrips=20, const bool bUnit=true)
std::vector< double > SubtractVec (std::vector< double > image0, std::vector< double > image1)
template<class T >
std::string ToString (const T &thing, int w=0, int p=0)
void ToTH2D (std::vector< double > vec, int nPlanes, int nStrips, std::string name, std::string title, TH2D *&Uview, TH2D *&Vview)
std::vector< double > ToVector (TH2D *hist, int nPlanes)
int Train (NtpStRecord *ntpst, std::vector< double > &avgNC, std::vector< double > &avgCCe, std::vector< double > &avgCCmu, std::vector< double > &avgCCtau, std::vector< double > &varNC, std::vector< double > &varCCe, std::vector< double > &varCCmu, std::vector< double > &varCCtau, std::vector< int > &numNC, std::vector< int > &numCCe, std::vector< int > &numCCmu, std::vector< int > &numCCtau, int &NCevts, int &eCCevts, int &muCCevts, int &tauCCevts, const int maxTrain=1000, const int nPlanes=20, const int nStrips=20, const int cutPlanes=40, const double eRecoMax=999., const double eTrueMax=999.)
void TrainPost (TFile *outFile, std::vector< double > &avgNC, std::vector< double > &avgCCe, std::vector< double > &avgCCmu, std::vector< double > &avgCCtau, std::vector< double > &varNC, std::vector< double > &varCCe, std::vector< double > &varCCmu, std::vector< double > &varCCtau, std::vector< int > &numNC, std::vector< int > &numCCe, std::vector< int > &numCCmu, std::vector< int > &numCCtau, int NCevts, int eCCevts, int muCCevts, int tauCCevts, const int nPlanes=20, const int nStrips=20)
void UnitVector (std::vector< double > &vec)
int VecIndex (int ipln, int istp, int nPlanes)
void WriteFile (TFile *outFile, std::vector< double > &avgNC, std::vector< double > &avgCCe, std::vector< double > &avgCCmu, std::vector< double > &avgCCtau, std::vector< double > &varNC, std::vector< double > &varCCe, std::vector< double > &varCCmu, std::vector< double > &varCCtau, std::vector< double > &pNChit, std::vector< double > &pCCehit, std::vector< double > &pCCmuhit, std::vector< double > &pCCtauhit, const int nPlanes, const int nStrips)

Enumeration Type Documentation

Enumerator:
vhsNC 
vhsCCe 
vhsCCmu 
vhsCCtau 
vhsUnknown 

Definition at line 28 of file VHS.h.

00028 { vhsNC=0, vhsCCe=1, vhsCCmu=2, vhsCCtau=3, vhsUnknown=4 };


Function Documentation

bool VHS::DrawEvent ( int  nEntry,
const char *  treeName,
int  nPlanes,
int  nStrips 
)

Definition at line 47 of file VHS.cxx.

References VHSevent::dCCe, VHSevent::dCCmu, VHSevent::dCCtau, VHSevent::dNC, FullVec(), GetEvtType(), VHSevent::GeV, Munits::GeV, VHSevent::ind, VHSevent::MC, VHSevent::MCenu, VHSevent::MCiaction, VHSevent::MCinu, VHSevent::NumPln, VHSevent::pCCe, VHSevent::pCCmu, VHSevent::pCCtau, VHSevent::pNC, VHSevent::Sigcor, VHSevent::stp, ToString(), ToTH2D(), vhsCCe, vhsCCmu, vhsCCtau, and vhsNC.

00049 {
00050   // Draw an event display for the given event from a TTree of VHSEvent
00051 
00052   TTree* theTree = (TTree*)gROOT->FindObject(treeName);
00053 
00054   if (!theTree) return false;
00055   const int kNumEntries = theTree->GetEntriesFast();
00056 
00057   if ( nEntry > kNumEntries ) return false;
00058 
00059   TCanvas* evtDisp = 0;
00060 
00061   if ( !gROOT->FindObject("VHSEvtDisp") ) {
00062     evtDisp = new TCanvas("VHSEvtDisp","VHS Event Display",0,0,800,600);
00063     evtDisp->Divide(2,2);
00064     evtDisp->GetPad(1)->SetPad(0.00,0.5,0.33,1.0);
00065     evtDisp->GetPad(2)->SetPad(0.33,0.5,1.00,1.0);
00066     evtDisp->GetPad(3)->SetPad(0.00,0.0,0.33,0.5);
00067     evtDisp->GetPad(4)->SetPad(0.33,0.0,1.00,0.5);
00068 
00069   } else {
00070     evtDisp = (TCanvas*)gROOT->FindObject("VHSEvtDisp");
00071     evtDisp->cd(1);
00072   }
00073   evtDisp->SetEditable(kTRUE);
00074 
00075   std::string nextCom = (nEntry+1 < kNumEntries)?
00076     std::string("VHS::DrawEvent(") + VHS::ToString(nEntry+1) + std::string(",\"") +
00077     std::string(theTree->GetName()) + "\"," + VHS::ToString(nPlanes) + std::string(",") +
00078     VHS::ToString(nStrips) + std::string(");")                                           :
00079     std::string("VHS::DrawEvent(") + VHS::ToString(kNumEntries-1) + std::string(",\"") +
00080     std::string(theTree->GetName()) + "\"," + VHS::ToString(nPlanes) + std::string(",") +
00081     VHS::ToString(nStrips) + std::string(");");
00082 
00083   std::string prevCom = (nEntry > 0)?
00084     std::string("VHS::DrawEvent(") + VHS::ToString(nEntry-1) + std::string(",\"") +
00085     std::string(theTree->GetName()) + "\"," + VHS::ToString(nPlanes) + std::string(",") +
00086     VHS::ToString(nStrips) + std::string(");")                                           :
00087     std::string("VHS::DrawEvent(0,\"") +
00088     std::string(theTree->GetName()) + "\"," + VHS::ToString(nPlanes) + std::string(",") +
00089     VHS::ToString(nStrips) + std::string(");");
00090 
00091   // Draw navigation buttons
00092   evtDisp->cd(3);
00093   evtDisp->GetPad(3)->GetListOfPrimitives()->Delete();
00094   evtDisp->GetPad(3)->Range(0,0,1,1);
00095 
00096   // Next button
00097   TButton* butNext =
00098     new TButton("Next Event",nextCom.c_str(),0.20,0.60,0.80,0.80);
00099   butNext->SetTextSize(0.3);
00100 
00101   // Previous button
00102   TButton* butPrev =
00103     new TButton("Prev Event",prevCom.c_str(),0.20,0.40,0.80,0.60);
00104   butPrev->SetTextSize(0.3);
00105     
00106   // Quit button
00107   TButton* butQuit = new TButton("Quit",".q",0.20,0.20,0.80,0.40);
00108   butQuit->SetTextSize(0.3);
00109 
00110   butNext->Draw();
00111   butPrev->Draw();
00112   butQuit->Draw();
00113 
00114   // Grab a copy of the event
00115   VHSevent* evt = new VHSevent();
00116   theTree->SetBranchAddress("VHSevent",&evt);
00117   theTree->GetEntry(nEntry);
00118 
00119   // Draw UZ view
00120   TH2D* uzImage = (TH2D*)gROOT->FindObject("UZImage");
00121   TH2D* vzImage = (TH2D*)gROOT->FindObject("VZImage");
00122 
00123   if ( uzImage ) { delete uzImage; uzImage = 0; }
00124   if ( vzImage ) { delete vzImage; vzImage = 0; }
00125 
00126   VHS::ToTH2D( VHS::FullVec( evt->stp, evt->ind, nPlanes, nStrips ),
00127                nPlanes, nStrips, "ZImage", "Z",
00128                uzImage, vzImage);
00129 
00130   evtDisp->cd(2);
00131   uzImage->SetXTitle("Plane");
00132   uzImage->SetYTitle("Strip");
00133   uzImage->SetStats(kFALSE);
00134   uzImage->Draw("COLZ");
00135 
00136   // Draw VZ view
00137   evtDisp->cd(4);
00138   vzImage->SetXTitle("Plane");
00139   vzImage->SetYTitle("Strip");
00140   vzImage->SetStats(kFALSE);
00141   vzImage->Draw("COLZ");
00142 
00143   // Draw Event information
00144   evtDisp->cd(1);
00145   evtDisp->GetPad(1)->GetListOfPrimitives()->Delete();
00146 
00147   bool   MC        = evt->MC;
00148   double MCenu     = evt->MCenu;
00149   int    MCiaction = evt->MCiaction;
00150   int    MCinu     = TMath::Abs(evt->MCinu);
00151   //int  MCrescode = evt->MCrescode;
00152 
00153   double dNC       = evt->dNC;
00154   double dCCe      = evt->dCCe;
00155   double dCCmu     = evt->dCCmu;
00156   double dCCtau    = evt->dCCtau;
00157 
00158   double pNC       = evt->pNC;
00159   double pCCe      = evt->pCCe;
00160   double pCCmu     = evt->pCCmu;
00161   double pCCtau    = evt->pCCtau;
00162 
00163   double GeV       = evt->GeV;
00164   double Sigcor    = evt->Sigcor;
00165 
00166   int    NumPln    = evt->NumPln;
00167 
00168   VHS::evtType evtt = VHS::GetEvtType( MCinu, MCiaction );
00169 
00170   // --> text boxes to dump info
00171   TPaveText* ptL = new TPaveText(0.02,0.02,0.5,0.45,"brNDC");
00172   ptL->SetLineWidth(2);
00173   ptL->SetBorderSize(1);
00174   ptL->SetTextSize(0.05);
00175   ptL->SetTextAlign(12);
00176 
00177   std::string li1L = std::string("dCCe = ")    + VHS::ToString(dCCe,7,5);
00178   std::string li2L = std::string("dCC#mu = ")  + VHS::ToString(dCCmu,7,5);
00179   std::string li3L = std::string("dCC#tau = ") + VHS::ToString(dCCtau,7,5);
00180   std::string li4L = std::string("dNC  = ")    + VHS::ToString(dNC,7,5);
00181   ptL->AddText(li1L.c_str());
00182   ptL->AddText(li2L.c_str());
00183   ptL->AddText(li3L.c_str());
00184   ptL->AddText(li4L.c_str());
00185   ptL->Draw();
00186 
00187   TPaveText* ptR = new TPaveText(0.5,0.02,0.98,0.45,"brNDC");
00188   ptR->SetLineWidth(2);
00189   ptR->SetBorderSize(1);
00190   ptR->SetTextSize(0.05);
00191   ptR->SetTextAlign(12);
00192 
00193   std::string li1R = std::string("pCCe = ")    + VHS::ToString(pCCe,7,5);
00194   std::string li2R = std::string("pCC#mu = ")  + VHS::ToString(pCCmu,7,5);
00195   std::string li3R = std::string("pCC#tau = ") + VHS::ToString(pCCtau,7,5);
00196   std::string li4R = std::string("pNC  = ")    + VHS::ToString(pNC,7,5);
00197   ptR->AddText(li1R.c_str());
00198   ptR->AddText(li2R.c_str());
00199   ptR->AddText(li3R.c_str());
00200   ptR->AddText(li4R.c_str());
00201   ptR->Draw();
00202 
00203   TPaveText* pt = new TPaveText(0.02,0.45,0.98,0.98,"brNDC");
00204   pt->SetLineWidth(2);
00205   pt->SetBorderSize(1);
00206   pt->SetTextSize(0.05);
00207   pt->SetTextAlign(12);
00208 
00209   std::string li1 =
00210     std::string("Event #") + VHS::ToString(nEntry+1) +
00211     std::string(" of ") + VHS::ToString(kNumEntries);
00212   std::string li2 = std::string("");
00213   if (evtt == vhsNC   ) li2 = "NC event ";
00214   if (evtt == vhsCCe  ) li2 = "e CCevent ";
00215   if (evtt == vhsCCmu ) li2 = "#mu CC event ";
00216   if (evtt == vhsCCtau) li2 = "#tau CC event ";
00217   li2 += (MC)?
00218     std::string("(") + VHS::ToString(NumPln) + std::string(" total planes)") :
00219     VHS::ToString(NumPln) + std::string(" total planes") ;
00220   std::string li3 = std::string("E_{vis} = ") + VHS::ToString(GeV,7,5) +
00221     std::string(" GeV (") + VHS::ToString(Sigcor,7,7) + std::string(" pe)");
00222   std::string li4 = (MC)?
00223     std::string("E_{#nu} = ") +
00224     VHS::ToString(MCenu,7,5)  +
00225     std::string(" GeV")       :
00226     "";
00227   std::string li5 = "";
00228   if ( dNC < dCCmu )
00229     li5 = std::string("Determined: NC event");
00230   if ( dCCmu < dNC )
00231     li5 = std::string("Determined: #mu CC event");
00232   //if ( dNC < dCCe && dNC < dCCmu && dNC < dCCtau )
00233   //li5 = std::string("Determined: NC event");
00234   //if ( dCCe < dNC && dCCe < dCCmu && dCCe < dCCtau )
00235   //li5 = std::string("Determined: e CC event");
00236   //if ( dCCmu < dNC && dCCmu < dCCe && dCCmu < dCCtau )
00237   //li5 = std::string("Determined: #mu CC event");
00238   //if ( dCCtau < dNC && dCCtau < dCCe && dCCtau < dCCmu )
00239   //li5 = std::string("Determined: #tau CC event");
00240   li5 += std::string(" (distance) ");
00241   pt->AddText(li1.c_str());
00242   pt->AddText(li2.c_str());
00243   pt->AddText(li3.c_str());
00244   if (MC) pt->AddText(li4.c_str());
00245   if ( (MC) ){
00246     if ((evtt == vhsNC    &&
00247          dNC > dCCmu       ) ||
00248         (evtt == vhsCCmu    &&
00249          dCCmu > dNC       )  )
00250          /*
00251          ( dNC    > dCCe  || dNC    > dCCmu || dNC    > dCCtau )) ||
00252         (evtt == vhsCCe   &&
00253          ( dCCe   > dNC   || dCCe   > dCCmu || dCCe   > dCCtau )) ||
00254         (evtt == vhsCCmu  &&
00255          ( dCCmu  > dNC   || dCCmu  > dCCe  || dCCmu  > dCCtau )) ||
00256         (evtt == vhsCCtau &&
00257          ( dCCtau > dNC   || dCCtau > dCCe  || dCCtau > dCCmu  ))  )
00258          */
00259       pt->AddText(li5.c_str())->SetTextColor(2);
00260     else
00261       pt->AddText(li5.c_str());
00262   } else {
00263     pt->AddText(li5.c_str());
00264   }
00265 
00266   std::string li6 = "";
00267   if ( pNC > pCCmu )
00268     li6 = std::string("Determined: NC event");
00269   if ( pCCmu > pNC )
00270     li6 = std::string("Determined: #mu CC event");
00271   /*
00272   if ( pNC > pCCe && pNC > pCCmu && pNC > pCCtau )
00273     li6 = std::string("Determined: NC event");
00274   if ( pCCe > pNC && pCCe > pCCmu && pCCe > pCCtau )
00275     li6 = std::string("Determined: e CC event");
00276   if ( pCCmu > pNC && pCCmu > pCCe && pCCmu > pCCtau )
00277     li6 = std::string("Determined: #mu CC event");
00278   if ( pCCtau > pNC && pCCtau > pCCe && pCCtau > pCCmu )
00279     li6 = std::string("Determined: #tau CC event");
00280   */
00281   li6 += std::string(" (loglikelihood) ");
00282   if ( (MC) ) {
00283     if ((evtt == vhsNC    &&
00284          pNC < pCCmu       ) ||
00285         (evtt == vhsCCmu  &&
00286          pCCmu < pNC       )  )
00287          /*
00288          ( pNC    < pCCe  || pNC    < pCCmu || pNC    < pCCtau )) ||
00289         (evtt == vhsCCe   &&
00290          ( pCCe   < pNC   || pCCe   < pCCmu || pCCe   < pCCtau )) ||
00291         (evtt == vhsCCmu  &&
00292          ( pCCmu  < pNC   || pCCmu  < pCCe  || pCCmu  < pCCtau )) ||
00293         (evtt == vhsCCtau &&
00294          ( pCCtau < pNC   || pCCtau < pCCe  || pCCtau < pCCmu  ))  )
00295          */
00296       pt->AddText(li6.c_str())->SetTextColor(2);
00297     else
00298       pt->AddText(li6.c_str());
00299   } else {
00300     pt->AddText(li6.c_str());
00301   }
00302 
00303   pt->Draw();
00304 
00305   delete evt; evt = 0;
00306 
00307   evtDisp->cd(3);
00308 
00309   evtDisp->SetEditable(kFALSE);
00310 
00311   return true;
00312 
00313 }

void VHS::FillDiscriminants ( NtpSREvent evt,
TClonesArray *  stp,
std::vector< double >  avgNC,
std::vector< double >  avgCCe,
std::vector< double >  avgCCmu,
std::vector< double >  avgCCtau,
std::vector< double >  varNC,
std::vector< double >  varCCe,
std::vector< double >  varCCmu,
std::vector< double >  varCCtau,
std::vector< double >  pNC,
std::vector< double >  pCCe,
std::vector< double >  pCCmu,
std::vector< double >  pCCtau,
const int  nPlanes,
const int  nStrips,
const bool  bUnit,
VHSevent *&  vhsevt 
)

Definition at line 317 of file VHS.cxx.

References NtpSRPlane::begu, NtpSRPlane::begv, VHSevent::dCCe, VHSevent::dCCmu, VHSevent::dCCtau, VHSevent::dNC, FullVec(), GetDistance(), GetImage(), GetLL(), VHSevent::ind, VHSevent::MinPlnU, VHSevent::MinPlnV, VHSevent::MOIThetaU, VHSevent::MOIThetaV, NtpSREvent::nstrip, VHSevent::pCCe, VHSevent::pCCmu, VHSevent::pCCtau, NtpSREvent::plane, VHSevent::pNC, NtpSREvent::stp, VHSevent::stp, and UnitVector().

Referenced by Skim().

00335 {
00336 
00337   // Gather our image and index sparse vectors
00338   std::vector<double> image;
00339   std::vector<int>    index;
00340   std::vector<double> theta;
00341   VHS::GetImage(evt->stp, evt->nstrip,
00342                 stp, nPlanes, nStrips, image, index, theta);
00343 
00344   vhsevt->stp = image;
00345   vhsevt->ind = index;
00346 
00347   vhsevt->MOIThetaU = theta[0];
00348   vhsevt->MOIThetaV = theta[1];
00349 
00350   // Get minimum plane in U & V views
00351   vhsevt->MinPlnU = evt->plane.begu;
00352   vhsevt->MinPlnV = evt->plane.begv;
00353 
00354   // Fill relative loglikelihoods
00355   std::vector<double> fullImg = FullVec(image,index,nPlanes,nStrips);
00356 
00357   vhsevt->pNC    =
00358     VHS::GetLL(fullImg, pNC   , avgNC   , varNC   );
00359   vhsevt->pCCe   =
00360     VHS::GetLL(fullImg, pCCe  , avgCCe  , varCCe  );
00361   vhsevt->pCCmu  =
00362     VHS::GetLL(fullImg, pCCmu , avgCCmu , varCCmu );
00363   vhsevt->pCCtau =
00364     VHS::GetLL(fullImg, pCCtau, avgCCtau, varCCtau);
00365 
00366   // Fill mean distances
00367   if (bUnit) {
00368     VHS::UnitVector(fullImg );
00369     VHS::UnitVector(avgNC   );
00370     VHS::UnitVector(avgCCe  );
00371     VHS::UnitVector(avgCCmu );
00372     VHS::UnitVector(avgCCtau);
00373   }
00374 
00375   vhsevt->dNC    = VHS::GetDistance( fullImg, avgNC    );
00376   vhsevt->dCCe   = VHS::GetDistance( fullImg, avgCCe   );
00377   vhsevt->dCCmu  = VHS::GetDistance( fullImg, avgCCmu  );
00378   vhsevt->dCCtau = VHS::GetDistance( fullImg, avgCCtau );
00379 
00380 }

std::vector< double > VHS::FindMedian ( std::vector< std::vector< double > >  pts,
std::vector< double >  avg,
bool  bUnit = true,
bool  bVerbose = true 
)

Definition at line 384 of file VHS.cxx.

References GetDistance(), OscPar::kDelta, and UnitVector().

00388 {
00389   const unsigned int nDim = avg.size(); // number of dimensions each point lives in
00390   const unsigned int nPts = pt.size();  // number of points for median-finding
00391   const double kEpsilon   = 1.e-6;      // precision to zero in calculations
00392   const double kDelta     = 1.e-4;      // cutoff for iterative median displacement
00393 
00394   if (bUnit) {
00395     VHS::UnitVector(avg);
00396     for (unsigned int i = 0; i < nPts; i++)
00397       VHS::UnitVector(pt[i]);
00398   }
00399 
00400   // Find the median via Weiszfeld iteration
00401   std::vector<double> median = avg;
00402   std::vector<double> lastMedian(nDim, 0.);
00403   double dDist    = 999.; // iterative difference in avg distance to iterative median
00404   double lastAvgD = 0.;   // average distance to the last iterative median
00405   double avgD     = 0.;   // average distance to the current iterative median
00406 
00407   // Find the average distance to our initial guess point (average)
00408   for (unsigned int i = 0; i < nPts; i++) {
00409     double dPt = VHS::GetDistance(median, pt[i]);
00410     avgD += dPt/nPts;
00411   }
00412   if (bVerbose)
00413     std::cout << "-- Iteration 0 : dItr = N/A : avgDist = " << avgD << " del(N/A)" << std::endl;
00414 
00415   // Iteratively move closer to the true median point
00416   double dItr     = 999.; // distance between last iterative median and new one
00417   int    nItr     = 0;    // number of iterations
00418   while( dItr > kDelta ) {
00419     lastMedian = median;
00420     for (unsigned int j = 0; j < nDim; j++) median[j] = 0.;
00421     lastAvgD = avgD;
00422     avgD = 0.;
00423 
00424     // Calculate our new median guess
00425     double norm = 0.;
00426     for (unsigned int i = 0; i < nPts; i++) {
00427       double dPt = VHS::GetDistance(lastMedian, pt[i]);
00428 
00429       if (!bUnit) norm += (dPt > kEpsilon)? 1./dPt : 1./kEpsilon;
00430       
00431       for (unsigned int j = 0; j < nDim; j++)
00432         median[j] += (dPt > kEpsilon)?
00433           pt[i][j]/dPt : pt[i][j]/TMath::Abs(kEpsilon) ;
00434     } // end for loop
00435     // Normalize our iterative median properly
00436     if (bUnit) VHS::UnitVector(median);
00437     else for (unsigned int j = 0; j < nDim; j++) median[j] *= 1./norm;
00438 
00439     dItr = VHS::GetDistance( lastMedian, median );
00440     // Find the average distance to our new guess point
00441     for (unsigned int i = 0; i < nPts; i++) {
00442       double dPt = VHS::GetDistance(median, pt[i]);
00443       avgD += dPt/nPts;
00444     }
00445 
00446     dDist = lastAvgD-avgD;
00447     nItr++;
00448     if (bVerbose)
00449       std::cout << "-- Iteration " << nItr << " : dItr = " << dItr
00450            << " : avgDist = " << avgD << " del(" << dDist
00451                 << ")" << std::endl;
00452   } // end while
00453   if (bVerbose) {
00454     double dMedAvg = VHS::GetDistance( median, avg );
00455     std::cout << "|median-average| = " << dMedAvg << std::endl;
00456   }
00457 
00458   return median;
00459 
00460 }

std::vector< double > VHS::FullVec ( std::vector< double >  image,
std::vector< int >  index,
int  nPlanes,
int  nStrips 
)

Definition at line 464 of file VHS.cxx.

Referenced by DrawEvent(), and FillDiscriminants().

00468 {
00469   std::vector<double> full( 2*nPlanes*nStrips, 0.);
00470 
00471   for (unsigned int i = 0; i < index.size(); i++)
00472     full[ index[i] ] = image[i];
00473 
00474   return full;
00475 
00476 }

double VHS::GetDistance ( std::vector< double >  vec0,
std::vector< double >  vec1 
)

Definition at line 480 of file VHS.cxx.

Referenced by FillDiscriminants(), and FindMedian().

00481 {
00482   // Return the Euclidean distance between the points described
00483   // by the input vectors
00484 
00485   if (vec0.size() != vec1.size()) return -1.;
00486 
00487   double dist = 0.;
00488 
00489   for (unsigned int i = 0; i < vec0.size(); i++)
00490     dist += TMath::Power( vec0[i]-vec1[i], 2. );
00491 
00492   return TMath::Sqrt( dist );
00493 
00494 }

VHS::evtType VHS::GetEvtType ( int  inu,
int  iaction 
)

Definition at line 498 of file VHS.cxx.

References vhsCCe, vhsCCmu, vhsCCtau, vhsNC, and vhsUnknown.

Referenced by DrawEvent(), Skim(), and Train().

00499 {
00500   if ( iaction == 0 ) return vhsNC; 
00501 
00502   if ( iaction == 1 && TMath::Abs(inu) == 12 ) return vhsCCe; 
00503 
00504   if ( iaction == 1 && TMath::Abs(inu) == 14 ) return vhsCCmu; 
00505 
00506   if ( iaction == 1 && TMath::Abs(inu) == 16 ) return vhsCCtau; 
00507 
00508   return vhsUnknown;
00509 
00510 }

void VHS::GetImage ( int *  index,
int  nstp,
TClonesArray *  stp,
int  nPlanes,
int  nStrips,
std::vector< double > &  image,
std::vector< int > &  vecInd,
std::vector< double > &  theta 
)

Definition at line 514 of file VHS.cxx.

References GetThetaAxis(), NtpSRStrip::ph0, NtpSRStrip::ph1, NtpSRStrip::plane, RotatePixel(), SeparateViews(), NtpSRPulseHeight::sigcor, NtpSRStrip::strip, and VecIndex().

Referenced by FillDiscriminants(), and Train().

00522 {
00523   // Return a standard vector of hit values representing
00524   // the event image for both U & V views
00525 
00526   // The images will be centered on strip ~nStrips/2 & plane ~nPlanes/4
00527 
00528   image.clear();    // empty the return vector of values
00529   vecInd.clear();   // empty the return vector of indices
00530   theta.clear();    // empty the return vector of event angle
00531   theta = std::vector<double>(2,-999.);
00532 
00533   // Separate the input index vector & stp array into u & v
00534   std::vector< std::vector<int> > uvindex(2); //uvindex index 0=u, 1=v
00535   VHS::SeparateViews(index,nstp,stp,uvindex[0],uvindex[1]);
00536 
00537   const int kOffsetInd = nPlanes*nStrips;
00538 
00539   for (unsigned int uv = 0; uv <= 1 ; uv++) {
00540     std::vector<int> ind = uvindex[uv];
00541 
00542     // Derive image rotation info
00543     std::vector<int> centerPix(2,0);
00544     VHS::GetThetaAxis(stp, uvindex[uv], theta[uv], centerPix);
00545 
00546     // Find the minimum plane in the view
00547     int minPln    = 999;
00548     for (unsigned int j = 0; j < ind.size(); j++) {
00549       NtpSRStrip* stpEntry = dynamic_cast<NtpSRStrip*>(stp->At(ind[j]));
00550       if (stpEntry) {
00551         int plane  = stpEntry->plane;
00552         int strip  = stpEntry->strip;
00553         VHS::RotatePixel(plane,strip,centerPix[0],centerPix[1],theta[uv]);
00554         if (plane < minPln) minPln = plane;
00555       } // end if stp ok
00556     } // end loop over strips
00557 
00558     // Find the end plane in our window
00559     int endPln = minPln + 2*nPlanes;
00560 
00561     // Find the weighted average plane in the view
00562     double stpXph    = 0.;
00563     double totSigcor = 0.;
00564     for (unsigned int j = 0; j < ind.size(); j++) {
00565       NtpSRStrip* stpEntry = dynamic_cast<NtpSRStrip*>(stp->At(ind[j]));
00566       if (stpEntry) {
00567         int plane  = stpEntry->plane;
00568         int strip  = stpEntry->strip;
00569         VHS::RotatePixel(plane,strip,centerPix[0],centerPix[1],theta[uv]);
00570         if (plane >= minPln && plane < endPln) {
00571           double sigcor = stpEntry->ph0.sigcor + stpEntry->ph1.sigcor;
00572 
00573           stpXph    += plane*sigcor;
00574           totSigcor += sigcor;
00575         } // end if strip is in our plane window
00576       } // end if stp ok
00577     } // end loop over strips
00578 
00579     // Find the weighted mean plane# in our window (~minos convention)
00580     // --> adjust the min & end plane numbers accordingly
00581     const int avgPln = TMath::FloorNint(stpXph/totSigcor);
00582     minPln = avgPln - TMath::FloorNint(0.25*2.*nPlanes);
00583     endPln = minPln + 2*nPlanes;
00584 
00585     // Find the weighted average strip in the view
00586     stpXph    = 0.;
00587     totSigcor = 0.;
00588     for (unsigned int j = 0; j < ind.size(); j++) {
00589       NtpSRStrip* stpEntry = dynamic_cast<NtpSRStrip*>(stp->At(ind[j]));
00590       if (stpEntry) {
00591         int plane  = stpEntry->plane;
00592         int strip  = stpEntry->strip;
00593         VHS::RotatePixel(plane,strip,centerPix[0],centerPix[1],theta[uv]);
00594         if (plane >= minPln && plane < endPln) {
00595           double sigcor = stpEntry->ph0.sigcor + stpEntry->ph1.sigcor;
00596 
00597           stpXph    += strip*sigcor;
00598           totSigcor += sigcor;
00599         } // end if strip is in our plane window
00600       } // end if stp ok
00601     } // end loop over strips
00602 
00603     // Find the weighted mean strip in our window
00604     const int avgStp = TMath::FloorNint(stpXph/totSigcor);
00605     const int minStp = avgStp - nStrips/2;
00606 
00607     // Loop over strips
00608     for (unsigned int i = 0; i < ind.size(); i++) {
00609       NtpSRStrip* stpEntry = dynamic_cast<NtpSRStrip*>(stp->At(ind[i]));
00610       if (stpEntry) {
00611         int plane  = stpEntry->plane;
00612         int strip  = stpEntry->strip;
00613         VHS::RotatePixel(plane,strip,centerPix[0],centerPix[1],theta[uv]);
00614         double sigcor = stpEntry->ph0.sigcor + stpEntry->ph1.sigcor;
00615 
00616         int iPln = (plane - minPln)/2;           // img coords
00617         int iStp = (strip - minStp);             // img coords
00618         int iVec = VecIndex(iPln,iStp,nPlanes);  // vec img coords
00619 
00620         if (iPln >= 0 && iPln < nPlanes &&
00621             iStp >= 0 && iStp < nStrips  ) {
00622           image.push_back(sigcor);
00623           if (uv == 0)
00624             vecInd.push_back(iVec);
00625           else
00626             vecInd.push_back(iVec+kOffsetInd);
00627         } // end if hit placement OK
00628       } // end if strip info OK
00629     } // end loop over strips
00630   } // end loop over u & v views
00631 
00632   return;
00633 
00634 }

double VHS::GetLL ( std::vector< double >  fullImage,
std::vector< double >  pHit,
std::vector< double >  avg,
std::vector< double >  var 
)

Definition at line 638 of file VHS.cxx.

Referenced by FillDiscriminants().

00642 {
00643   // Given a set of averages and variances,
00644   // calculate the Gaussian probability of the collection of
00645   // particular hit values convoluted with the probability of the
00646   // particular hit pattern
00647 
00648   double ll = 0.;
00649   const double kCutoff   = 1.e-10;
00650   const double kLnCutoff = TMath::Log10(kCutoff);
00651 
00652   for (unsigned int i = 0; i < fullImage.size(); i++) {
00653     double x     = fullImage[i];
00654     double mean  = avg[i];
00655     double sigma = TMath::Sqrt( var[i] );
00656 
00657     if ( x < kCutoff || pHit[i] < kCutoff ) 
00658       ll += TMath::Log10(1.-pHit[i]);
00659     else {
00660       ll += TMath::Log10(pHit[i]);
00661       ll += (TMath::Gaus(x,mean,sigma,true) > kCutoff)?
00662         TMath::Log10(TMath::Gaus(x,mean,sigma,true)) : kLnCutoff;
00663       ll -= (TMath::Gaus(mean,mean,sigma,true) > kCutoff)?
00664         TMath::Log10(TMath::Gaus(mean,mean,sigma,true)) : kLnCutoff;
00665     } // end if values within ranges
00666   } // end loop over image entries
00667 
00668   return ll;
00669 
00670 }

int VHS::GetPlane ( int  vecInd,
int  nPlanes 
)

Definition at line 674 of file VHS.cxx.

Referenced by TrackSegmentCamAtNu::AddCluster(), ShowerSegmentCamAtNu::AddCluster(), TrackSegmentCam::AddCluster(), AlgTrackSRList::AddClustersToTracks(), TrackCamAtNu::AddHit(), TrackCam::AddHit(), UberDST::Ana(), Anp::FillShortEvent::BVD_BestLine(), AlgTrackSRList::CheckForBadClusters(), Managed::HitManager::ClearXTalk(), Anp::PlotPmt::Collect(), NtpSRModule::FillNtpTrackTime(), UberModule::FillNtpTrackTime(), UberModuleLite::FillNtpTrackTime(), Managed::HitManager::FindHit(), AlgShowerCam::FindShowerVertex(), AlgAtmosShowerList::FormDummyTracks(), Anp::FillShortVar::Get(), ShowerSegmentCamAtNu::GetBegDir(), TrackSegmentCam::GetBegDir(), TrackSegmentCamAtNu::GetBegDir(), ShowerSegmentCamAtNu::GetBegStrip(), TrackSegmentCam::GetBegTPos(), TrackSegmentCamAtNu::GetBegTPos(), ShowerSegmentCamAtNu::GetBegTPos(), ObjShowerAtNu::GetDir(), TrackSegmentCam::GetEndDir(), TrackSegmentCamAtNu::GetEndDir(), ShowerSegmentCamAtNu::GetEndDir(), ShowerSegmentCamAtNu::GetEndStrip(), TrackSegmentCamAtNu::GetEndTPos(), ShowerSegmentCamAtNu::GetEndTPos(), TrackSegmentCam::GetEndTPos(), CandSliceHandle::GetNPlane(), Anp::FillShortEvent::GetPlaneStrips(), LISummaryModule::GetSummaryBlocks(), Anp::FillShortVar::GetTrackMap(), NtpSRBleachFiller::initStraightPHFraction(), AlgTrackSRList::MergeTracks(), Reco::ByPlane::operator()(), VtxClusterList::Process(), UberModule::Reco(), UberModuleLite::Reco(), VtxClusterList::Report(), Anp::FillSnarl::Run(), AlgAtNuRecoMCTruth::RunAlg(), AlgShowerAtNu::RunAlg(), AlgFarDetEvent::RunAlg(), AlgShowerSR::SetUV(), AlgTrack::SetUVZ(), AlgTrackSRList::SpectrometerTracking(), AlgClusterSRList::StripKeyFromPlane(), AlgSubShowerSRList::StripKeyFromPlane(), StripSRKeyFromPlane(), AlgSubShowerSRList::TestOverLap(), and ToTH2D().

00675 {
00676   // Returns relative plane # (*not* MINOS global plane convention)
00677   return (nPlanes != 0)? vecInd % nPlanes : -1;
00678 }

int VHS::GetStrip ( int  vecInd,
int  nPlanes 
)
void VHS::GetThetaAxis ( TClonesArray *  stp,
std::vector< int >  index,
double &  theta,
std::vector< int > &  center 
)

Definition at line 690 of file VHS.cxx.

References NtpSRStrip::ph0, NtpSRStrip::ph1, NtpSRStrip::plane, NtpSRPulseHeight::sigcor, and NtpSRStrip::strip.

Referenced by GetImage().

00694 {
00695   const unsigned int kNhits = index.size();
00696 
00697   // Find the weighted average plane & strip in the view
00698   // And the minimum plane/strip
00699   double stpXph    = 0.;
00700   double plnXph    = 0.;
00701   double totSigcor = 0.;
00702   int    minPln    = 999;
00703   int    minPlnStp = 999;
00704   for (unsigned int j = 0; j < kNhits; j++) {
00705     NtpSRStrip* stpEntry = dynamic_cast<NtpSRStrip*>(stp->At(index[j]));
00706     if (stpEntry) {
00707       int    plane  = stpEntry->plane;
00708       int    strip  = stpEntry->strip;
00709       double sigcor = stpEntry->ph0.sigcor + stpEntry->ph1.sigcor;
00710 
00711       plnXph    += plane*sigcor;
00712       stpXph    += strip*sigcor;
00713       totSigcor += sigcor;
00714 
00715       if (plane < minPln) {
00716         minPln    = plane;
00717         minPlnStp = strip;
00718       } // end if minimum plane number
00719     } // end if stp ok
00720   } // end loop over strips
00721 
00722   center[0] = (int)TMath::Floor(plnXph/totSigcor);
00723   center[1] = (int)TMath::Floor(stpXph/totSigcor);
00724 
00725   double relPln = plnXph/totSigcor - 1.*minPln;
00726   double relStp = stpXph/totSigcor - 1.*minPlnStp;
00727 
00728   theta = (TMath::Abs(relPln) > 1.e-10)?
00729     TMath::ATan( relStp/relPln ) : 0.;
00730 
00731   return;
00732 }

std::vector< std::string > VHS::ParseNmList ( const char *  cstr  ) 

Definition at line 736 of file VHS.cxx.

References len.

00737 {
00738   // Turn the input space delimited cstring into a vector of strings
00739 
00740   std::string strList = cstr;
00741 
00742   std::vector<std::string> strVec;
00743 
00744   int len = strList.length();
00745 
00746   unsigned int lastLoc = 0;
00747   unsigned int loc = strList.find(" ",lastLoc);
00748   while (lastLoc != std::string::npos) {
00749     int subLen = (loc != std::string::npos)? loc-lastLoc : len-lastLoc;
00750     std::string sub = strList.substr(lastLoc, subLen);
00751     if (sub.length() > 0) strVec.push_back(sub);
00752     lastLoc = (loc != std::string::npos)? loc+1 : loc;
00753     loc = strList.find(" ",lastLoc);
00754   }
00755 
00756   return strVec;
00757 }

void VHS::ReadFile ( TFile *  inFile,
std::vector< double > &  avgNC,
std::vector< double > &  avgCCe,
std::vector< double > &  avgCCmu,
std::vector< double > &  avgCCtau,
std::vector< double > &  varNC,
std::vector< double > &  varCCe,
std::vector< double > &  varCCmu,
std::vector< double > &  varCCtau,
std::vector< double > &  pNChit,
std::vector< double > &  pCCehit,
std::vector< double > &  pCCmuhit,
std::vector< double > &  pCCtauhit,
const int  nPlanes 
)

Definition at line 761 of file VHS.cxx.

References ToVector().

00775 {
00776   //
00777   // inFile must be open and valid
00778   //
00779   // Translate averages & completed variances from TH2D and write
00780   // to vectors for use
00781   //
00782 
00783   // Prepare to write our training info to output file
00784   inFile->cd();
00785 
00786   // NC hit prob hist
00787   const std::string UpNCName  = "UpNCHit";
00788   const std::string VpNCName  = "VpNCHit";
00789   TH2D* UpNCHist = (TH2D*)inFile->Get(UpNCName.c_str());
00790   TH2D* VpNCHist = (TH2D*)inFile->Get(VpNCName.c_str());
00791   std::vector<double> UpNChit = VHS::ToVector( UpNCHist, nPlanes);
00792   std::vector<double> VpNChit = VHS::ToVector( VpNCHist, nPlanes);
00793   pNChit.clear();
00794   for (unsigned int i = 0; i < UpNChit.size(); i++)
00795     pNChit.push_back(UpNChit[i]);
00796   for (unsigned int i = 0; i < VpNChit.size(); i++)
00797     pNChit.push_back(VpNChit[i]);
00798 
00799   // CCe hit prob hist
00800   const std::string UpCCeName  = "UpCCeHit";
00801   const std::string VpCCeName  = "VpCCeHit";
00802   TH2D* UpCCeHist = (TH2D*)inFile->Get(UpCCeName.c_str());
00803   TH2D* VpCCeHist = (TH2D*)inFile->Get(VpCCeName.c_str());
00804   std::vector<double> UpCCehit = VHS::ToVector( UpCCeHist, nPlanes);
00805   std::vector<double> VpCCehit = VHS::ToVector( VpCCeHist, nPlanes);
00806   pCCehit.clear();
00807   for (unsigned int i = 0; i < UpCCehit.size(); i++)
00808     pCCehit.push_back(UpCCehit[i]);
00809   for (unsigned int i = 0; i < VpCCehit.size(); i++)
00810     pCCehit.push_back(VpCCehit[i]);
00811 
00812   // CCmu hit prob hist
00813   const std::string UpCCmuName  = "UpCCmuHit";
00814   const std::string VpCCmuName  = "VpCCmuHit";
00815   TH2D* UpCCmuHist = (TH2D*)inFile->Get(UpCCmuName.c_str());
00816   TH2D* VpCCmuHist = (TH2D*)inFile->Get(VpCCmuName.c_str());
00817   std::vector<double> UpCCmuhit = VHS::ToVector( UpCCmuHist, nPlanes);
00818   std::vector<double> VpCCmuhit = VHS::ToVector( VpCCmuHist, nPlanes);
00819   pCCmuhit.clear();
00820   for (unsigned int i = 0; i < UpCCmuhit.size(); i++)
00821     pCCmuhit.push_back(UpCCmuhit[i]);
00822   for (unsigned int i = 0; i < VpCCmuhit.size(); i++)
00823     pCCmuhit.push_back(VpCCmuhit[i]);
00824 
00825   // CCtau hit prob hist
00826   const std::string UpCCtauName  = "UpCCtauHit";
00827   const std::string VpCCtauName  = "VpCCtauHit";
00828   TH2D* UpCCtauHist = (TH2D*)inFile->Get(UpCCtauName.c_str());
00829   TH2D* VpCCtauHist = (TH2D*)inFile->Get(VpCCtauName.c_str());
00830   std::vector<double> UpCCtauhit = VHS::ToVector( UpCCtauHist, nPlanes);
00831   std::vector<double> VpCCtauhit = VHS::ToVector( VpCCtauHist, nPlanes);
00832   pCCtauhit.clear();
00833   for (unsigned int i = 0; i < UpCCtauhit.size(); i++)
00834     pCCtauhit.push_back(UpCCtauhit[i]);
00835   for (unsigned int i = 0; i < VpCCtauhit.size(); i++)
00836     pCCtauhit.push_back(VpCCtauhit[i]);
00837 
00838   // NC Avg hist
00839   const std::string UavgNCName  = "UAvgNCImg";
00840   const std::string VavgNCName  = "VAvgNCImg";
00841   TH2D* UavgNCHist = (TH2D*)inFile->Get(UavgNCName.c_str());
00842   TH2D* VavgNCHist = (TH2D*)inFile->Get(VavgNCName.c_str());
00843   std::vector<double> UavgNC = VHS::ToVector( UavgNCHist, nPlanes);
00844   std::vector<double> VavgNC = VHS::ToVector( VavgNCHist, nPlanes);
00845   avgNC.clear();
00846   for (unsigned int i = 0; i < UavgNC.size(); i++)
00847     avgNC.push_back(UavgNC[i]);
00848   for (unsigned int i = 0; i < VavgNC.size(); i++)
00849     avgNC.push_back(VavgNC[i]);
00850 
00851   // CCe Avg hist
00852   const std::string UavgCCeName  = "UAvgCCeImg";
00853   const std::string VavgCCeName  = "VAvgCCeImg";
00854   TH2D* UavgCCeHist = (TH2D*)inFile->Get(UavgCCeName.c_str());
00855   TH2D* VavgCCeHist = (TH2D*)inFile->Get(VavgCCeName.c_str());
00856   std::vector<double> UavgCCe = VHS::ToVector( UavgCCeHist, nPlanes);
00857   std::vector<double> VavgCCe = VHS::ToVector( VavgCCeHist, nPlanes);
00858   avgCCe.clear();
00859   for (unsigned int i = 0; i < UavgCCe.size(); i++)
00860     avgCCe.push_back(UavgCCe[i]);
00861   for (unsigned int i = 0; i < VavgCCe.size(); i++)
00862     avgCCe.push_back(VavgCCe[i]);
00863 
00864   // CCmu Avg hist
00865   const std::string UavgCCmuName  = "UAvgCCmuImg";
00866   const std::string VavgCCmuName  = "VAvgCCmuImg";
00867   TH2D* UavgCCmuHist = (TH2D*)inFile->Get(UavgCCmuName.c_str());
00868   TH2D* VavgCCmuHist = (TH2D*)inFile->Get(VavgCCmuName.c_str());
00869   std::vector<double> UavgCCmu = VHS::ToVector( UavgCCmuHist, nPlanes);
00870   std::vector<double> VavgCCmu = VHS::ToVector( VavgCCmuHist, nPlanes);
00871   avgCCmu.clear();
00872   for (unsigned int i = 0; i < UavgCCmu.size(); i++)
00873     avgCCmu.push_back(UavgCCmu[i]);
00874   for (unsigned int i = 0; i < VavgCCmu.size(); i++)
00875     avgCCmu.push_back(VavgCCmu[i]);
00876 
00877   // CCtau Avg hist
00878   const std::string UavgCCtauName  = "UAvgCCtauImg";
00879   const std::string VavgCCtauName  = "VAvgCCtauImg";
00880   TH2D* UavgCCtauHist = (TH2D*)inFile->Get(UavgCCtauName.c_str());
00881   TH2D* VavgCCtauHist = (TH2D*)inFile->Get(VavgCCtauName.c_str());
00882   std::vector<double> UavgCCtau = VHS::ToVector( UavgCCtauHist, nPlanes);
00883   std::vector<double> VavgCCtau = VHS::ToVector( VavgCCtauHist, nPlanes);
00884   avgCCtau.clear();
00885   for (unsigned int i = 0; i < UavgCCtau.size(); i++)
00886     avgCCtau.push_back(UavgCCtau[i]);
00887   for (unsigned int i = 0; i < VavgCCtau.size(); i++)
00888     avgCCtau.push_back(VavgCCtau[i]);
00889 
00890   // NC variance hist
00891   const std::string UNCVarName  = "UNCVar";
00892   const std::string VNCVarName  = "VNCVar";
00893   TH2D* UNCVarHist = (TH2D*)inFile->Get(UNCVarName.c_str());
00894   TH2D* VNCVarHist = (TH2D*)inFile->Get(VNCVarName.c_str());
00895   std::vector<double> UNCVar = VHS::ToVector( UNCVarHist, nPlanes);
00896   std::vector<double> VNCVar = VHS::ToVector( VNCVarHist, nPlanes);
00897   varNC.clear();
00898   for (unsigned int i = 0; i < UNCVar.size(); i++)
00899     varNC.push_back(UNCVar[i]);
00900   for (unsigned int i = 0; i < VpNChit.size(); i++)
00901     varNC.push_back(VNCVar[i]);
00902 
00903   // CCe variance hist
00904   const std::string UCCeVarName  = "UCCeVar";
00905   const std::string VCCeVarName  = "VCCeVar";
00906   TH2D* UCCeVarHist = (TH2D*)inFile->Get(UCCeVarName.c_str());
00907   TH2D* VCCeVarHist = (TH2D*)inFile->Get(VCCeVarName.c_str());
00908   std::vector<double> UCCeVar = VHS::ToVector( UCCeVarHist, nPlanes);
00909   std::vector<double> VCCeVar = VHS::ToVector( VCCeVarHist, nPlanes);
00910   varCCe.clear();
00911   for (unsigned int i = 0; i < UCCeVar.size(); i++)
00912     varCCe.push_back(UCCeVar[i]);
00913   for (unsigned int i = 0; i < VpCCehit.size(); i++)
00914     varCCe.push_back(VCCeVar[i]);
00915 
00916   // CCmu variance hist
00917   const std::string UCCmuVarName  = "UCCmuVar";
00918   const std::string VCCmuVarName  = "VCCmuVar";
00919   TH2D* UCCmuVarHist = (TH2D*)inFile->Get(UCCmuVarName.c_str());
00920   TH2D* VCCmuVarHist = (TH2D*)inFile->Get(VCCmuVarName.c_str());
00921   std::vector<double> UCCmuVar = VHS::ToVector( UCCmuVarHist, nPlanes);
00922   std::vector<double> VCCmuVar = VHS::ToVector( VCCmuVarHist, nPlanes);
00923   varCCmu.clear();
00924   for (unsigned int i = 0; i < UCCmuVar.size(); i++)
00925     varCCmu.push_back(UCCmuVar[i]);
00926   for (unsigned int i = 0; i < VpCCmuhit.size(); i++)
00927     varCCmu.push_back(VCCmuVar[i]);
00928 
00929   // CCtau variance hist
00930   const std::string UCCtauVarName  = "UCCtauVar";
00931   const std::string VCCtauVarName  = "VCCtauVar";
00932   TH2D* UCCtauVarHist = (TH2D*)inFile->Get(UCCtauVarName.c_str());
00933   TH2D* VCCtauVarHist = (TH2D*)inFile->Get(VCCtauVarName.c_str());
00934   std::vector<double> UCCtauVar = VHS::ToVector( UCCtauVarHist, nPlanes);
00935   std::vector<double> VCCtauVar = VHS::ToVector( VCCtauVarHist, nPlanes);
00936   varCCtau.clear();
00937   for (unsigned int i = 0; i < UCCtauVar.size(); i++)
00938     varCCtau.push_back(UCCtauVar[i]);
00939   for (unsigned int i = 0; i < VpCCtauhit.size(); i++)
00940     varCCtau.push_back(VCCtauVar[i]);
00941 
00942   inFile->Close();
00943 
00944   return;
00945 
00946 }

void VHS::RotatePixel ( int &  plane,
int &  strip,
const int  avgPlane,
const int  avgStrip,
const double  theta 
)

Definition at line 950 of file VHS.cxx.

References MuELoss::e.

Referenced by GetImage().

00953 {
00954   // rotate around pixel center to new coords (pln, stp);
00955   // --> round floats to nearest integer
00956 
00957   if ( TMath::Abs(theta) < 1.e-10 ) return; // ~zero rotation
00958 
00959   const double cosTheta = TMath::Cos(-theta); // rotate opposite dir as given
00960   const double sinTheta = TMath::Sin(-theta); // rotate opposite dir as given
00961 
00962   double newPlnD = cosTheta*(plane-avgPlane) + sinTheta*(strip-avgStrip);
00963   double newStpD = cosTheta*(strip-avgStrip) - sinTheta*(plane-avgPlane);
00964 
00965   if (newStpD-TMath::Floor(newStpD) < 0.5) newStpD = TMath::Floor(newStpD);
00966   else newStpD = TMath::Ceil(newStpD);
00967 
00968   if (newPlnD-TMath::Floor(newPlnD) < 0.5) newPlnD = TMath::Floor(newPlnD);
00969   else newPlnD = TMath::Ceil(newPlnD);
00970 
00971   strip = avgStrip + (int)newStpD;
00972   plane = avgPlane + (int)newPlnD;
00973 
00974   return;
00975 
00976 }

void VHS::SeparateViews ( int *  index,
int  nstp,
TClonesArray *  stp,
std::vector< int > &  ustp,
std::vector< int > &  vstp 
)

Definition at line 980 of file VHS.cxx.

References PlaneView::kU, PlaneView::kV, and NtpSRStrip::planeview.

Referenced by GetImage().

00982 {
00983   // Sort hit strips into separate views
00984 
00985   for (int i = 0; i < nstp; i++) {
00986     int ind = index[i];
00987     NtpSRStrip* strip = dynamic_cast<NtpSRStrip*>(stp->At(ind));
00988     if (strip) {
00989       unsigned int pv = strip->planeview;
00990       if (pv == PlaneView::kU) ustp.push_back(ind);
00991       if (pv == PlaneView::kV) vstp.push_back(ind);
00992     }
00993   }
00994   return;
00995 }

std::vector< VHSevent * > VHS::Skim ( NtpStRecord ntpst,
std::vector< double >  avgNC,
std::vector< double >  avgCCe,
std::vector< double >  avgCCmu,
std::vector< double >  avgCCtau,
std::vector< double >  varNC,
std::vector< double >  varCCe,
std::vector< double >  varCCmu,
std::vector< double >  varCCtau,
std::vector< double >  pNC,
std::vector< double >  pCCe,
std::vector< double >  pCCmu,
std::vector< double >  pCCtau,
const int  nPlanes = 20,
const int  nStrips = 20,
const bool  bUnit = true 
)

Definition at line 999 of file VHS.cxx.

References NtpSRPlane::beg, NtpSRPlane::end, NtpStRecord::evt, NtpStRecord::evthdr, FillDiscriminants(), GetEvtType(), RecRecordImp< T >::GetHeader(), NtpSRStripPulseHeight::gev, VHSevent::GeV, NtpMCTruth::iaction, NtpMCTruth::inu, NtpMCTruth::iresonance, VHSevent::MC, NtpStRecord::mc, VHSevent::MCenu, VHSevent::MCiaction, VHSevent::MCinu, VHSevent::MCQsq, VHSevent::MCrescode, VHSevent::MCWsq, VHSevent::MCx, VHSevent::MCy, NtpTHEvent::neumc, NtpSREventSummary::nevent, NtpSREvent::nstrip, VHSevent::NumPln, VHSevent::NumStp, NtpMCTruth::p4neu, NtpSREvent::ph, NtpSREvent::plane, NtpMCTruth::q2, NtpSRPulseHeight::sigcor, VHSevent::Sigcor, NtpStRecord::stp, NtpStRecord::thevt, vhsUnknown, NtpMCTruth::w2, NtpMCTruth::x, and NtpMCTruth::y.

01015 {
01016   // Inputs: NtpStRecord, average images, variances,
01017   //         vector of raw hit probabilities, window size (pln x stp),
01018   //         flag signalling whether to compute distances on unit or
01019   //         absolute image vectors
01020   // Output: Vector of pointers to complete VHSevent records
01021   //
01022 
01023   // Define the vector we will return
01024   std::vector< VHSevent* > vhslist;
01025 
01026   // If we have no record, do nothing
01027   if ( !ntpst ) return vhslist;
01028 
01029   const NtpSREventSummary& evthdr = ntpst->evthdr;
01030   TClonesArray*            evt    = ntpst->evt;
01031   TClonesArray*            stp    = ntpst->stp;
01032   TClonesArray*            mc     = ntpst->mc;
01033   TClonesArray*            thevt  = ntpst->thevt;
01034 
01035   // Event quantities
01036   if (evt && stp) {
01037     const unsigned int nEvt = evthdr.nevent; // evt->GetEntries();
01038 
01039     for (unsigned int i = 0; i < nEvt; i++){
01040       NtpSREvent* evtEntry = dynamic_cast<NtpSREvent*>( evt->At(i) );
01041       if ( !evtEntry ) continue;
01042 
01043       // Define our VHSevent to fill
01044       VHSevent* vhsevt = new VHSevent( ntpst->GetHeader() ) ;
01045       vhsevt->SetName("VHSevent");
01046       vhsevt->SetTitle("VHS Event");
01047 
01048       // Fill NtpSREvent info
01049       vhsevt->GeV    = evtEntry->ph.gev;
01050       vhsevt->NumPln = evtEntry->plane.end - evtEntry->plane.beg;
01051       vhsevt->NumStp = evtEntry->nstrip;
01052       vhsevt->Sigcor = evtEntry->ph.sigcor;
01053 
01054       VHS::evtType evtt = vhsUnknown;
01055 
01056       // Fill MC information
01057       if ( mc && thevt ) {
01058         NtpTHEvent*  thevtEntry  = dynamic_cast<NtpTHEvent*>(thevt->At(i));
01059         NtpMCTruth*  mcEntry     = (thevtEntry)?
01060           dynamic_cast<NtpMCTruth*>( mc->At( (thevtEntry->neumc) ) ) : 0;
01061         if (mcEntry) {
01062           vhsevt->MCenu     = mcEntry->p4neu[3];
01063           vhsevt->MCinu     = mcEntry->inu;
01064           vhsevt->MCiaction = mcEntry->iaction;
01065           vhsevt->MCrescode = mcEntry->iresonance;
01066           vhsevt->MCQsq     = mcEntry->q2;
01067           vhsevt->MCWsq     = mcEntry->w2;
01068           vhsevt->MCx       = mcEntry->x;
01069           vhsevt->MCy       = mcEntry->y;
01070 
01071           evtt = VHS::GetEvtType(mcEntry->inu, mcEntry->iaction);
01072         } // end if mcEntry
01073       } // end if mc && thevt
01074       vhsevt->MC = (evtt != vhsUnknown);
01075 
01076       VHS::FillDiscriminants(evtEntry, stp,
01077                              avgNC, avgCCe, avgCCmu, avgCCtau,
01078                              varNC, varCCe, varCCmu, varCCtau,
01079                              pNC,   pCCe,   pCCmu,   pCCtau,
01080                              nPlanes, nStrips, bUnit, vhsevt );
01081 
01082       vhslist.push_back( vhsevt ); // add this event to the list
01083 
01084     } // end loop over event entries
01085   } // end if evt and stp arrays exist
01086 
01087   return vhslist;
01088 }

std::vector< double > VHS::SubtractVec ( std::vector< double >  image0,
std::vector< double >  image1 
)

Definition at line 1092 of file VHS.cxx.

01094 {
01095   // 'nuff said
01096 
01097   std::vector<double> subt(image0.size(), 0.);
01098 
01099   if (image0.size() != image1.size() ) return subt;
01100 
01101   for (unsigned int i = 0; i < image0.size(); i++)
01102     subt[i] = image0[i] - image1[i];
01103 
01104   return subt;
01105 
01106 }

template<class T >
std::string VHS::ToString ( const T &  thing,
int  w = 0,
int  p = 0 
) [inline]

Definition at line 1111 of file VHS.cxx.

Referenced by FiltTriggerPrescale::Config(), FiltTriggerPrescale::DefaultConfig(), and DrawEvent().

01112 {
01113 
01114   std::ostringstream os;
01115   os << std::setw(w) << std::setprecision(p) << thing;
01116   return os.str();
01117 
01118 }

void VHS::ToTH2D ( std::vector< double >  vec,
int  nPlanes,
int  nStrips,
std::string  name,
std::string  title,
TH2D *&  Uview,
TH2D *&  Vview 
)

Definition at line 1122 of file VHS.cxx.

References GetPlane(), and GetStrip().

Referenced by DrawEvent(), and WriteFile().

01126 {
01127   // 'nuff said
01128 
01129   std::string Uname  = "U" + name;
01130   std::string Utitle = "U" + title;
01131 
01132   std::string Vname  = "V" + name;
01133   std::string Vtitle = "V" + title;
01134 
01135   //if (Uview) delete Uview;
01136   //if (Vview) delete Vview;
01137 
01138   Uview = new TH2D(Uname.c_str(), Utitle.c_str(),
01139                    nPlanes, 0, nPlanes,
01140                    nStrips, 0, nStrips);
01141   Vview = new TH2D(Vname.c_str(), Vtitle.c_str(),
01142                    nPlanes, 0, nPlanes,
01143                    nStrips, 0, nStrips);
01144   const unsigned int kOffset = nPlanes*nStrips;
01145   for (unsigned int i = 0; i < vec.size(); i++){
01146     if ( i < kOffset )
01147       Uview->SetBinContent( VHS::GetPlane(i,nPlanes)+1,
01148                             VHS::GetStrip(i,nPlanes)+1,
01149                             vec[i]                  );
01150     else
01151       Vview->SetBinContent( VHS::GetPlane(i-kOffset,nPlanes)+1,
01152                             VHS::GetStrip(i-kOffset,nPlanes)+1,
01153                             vec[i]                  );
01154   } // end loop over vec entries
01155 
01156   return;
01157 
01158 }

std::vector< double > VHS::ToVector ( TH2D *  hist,
int  nPlanes 
)

Definition at line 1162 of file VHS.cxx.

References VecIndex().

Referenced by ReadFile().

01163 {
01164   // 'nuff said
01165 
01166   const int nPln  = hist->GetNbinsX();
01167   const int nStp  = hist->GetNbinsY();
01168   const int nBins = nPln*nStp;
01169 
01170   std::vector<double> vec(nBins,0.);
01171 
01172   if (nPln != nPlanes) return vec;
01173 
01174   for (int ipln = 0; ipln < nPln; ipln++)
01175     for (int istp = 0; istp < nStp; istp++)
01176       vec[VecIndex(ipln,istp,nPlanes)] =
01177         hist->GetBinContent(ipln+1, istp+1);
01178 
01179   return vec;
01180 
01181 }

int VHS::Train ( NtpStRecord ntpst,
std::vector< double > &  avgNC,
std::vector< double > &  avgCCe,
std::vector< double > &  avgCCmu,
std::vector< double > &  avgCCtau,
std::vector< double > &  varNC,
std::vector< double > &  varCCe,
std::vector< double > &  varCCmu,
std::vector< double > &  varCCtau,
std::vector< int > &  numNC,
std::vector< int > &  numCCe,
std::vector< int > &  numCCmu,
std::vector< int > &  numCCtau,
int &  NCevts,
int &  eCCevts,
int &  muCCevts,
int &  tauCCevts,
const int  maxTrain = 1000,
const int  nPlanes = 20,
const int  nStrips = 20,
const int  cutPlanes = 40,
const double  eRecoMax = 999.,
const double  eTrueMax = 999. 
)

Definition at line 1185 of file VHS.cxx.

References NtpSRPlane::beg, NtpSRPlane::end, NtpStRecord::evt, NtpStRecord::evthdr, GetEvtType(), GetImage(), NtpSRStripPulseHeight::gev, NtpMCTruth::iaction, NtpMCTruth::inu, NtpStRecord::mc, NtpTHEvent::neumc, NtpSREventSummary::nevent, NtpSREvent::nstrip, NtpMCTruth::p4neu, NtpSREvent::ph, NtpSREvent::plane, NtpSREvent::stp, NtpStRecord::stp, NtpStRecord::thevt, vhsCCe, vhsCCmu, vhsCCtau, vhsNC, and vhsUnknown.

01208 {
01209   // Input : NtpStRecord, average images to modify, variances to
01210   //         modify, vectors of number of hits to each pixel to
01211   //         modify, #planes in each view window, #strips in each view
01212   //         window
01213   // Output: Modified average images & modified covariance matrices
01214   //
01215 
01216   int nEvtUsed = 0;
01217 
01218   // If we have no record, do nothing
01219   if ( !ntpst ) return nEvtUsed;
01220 
01221   const NtpSREventSummary& evthdr = ntpst->evthdr;
01222   TClonesArray*            evt    = ntpst->evt;
01223   TClonesArray*            stp    = ntpst->stp;
01224   TClonesArray*            mc     = ntpst->mc;
01225   TClonesArray*            thevt  = ntpst->thevt;
01226 
01227   // Event quantities
01228   if (evt && stp && mc && thevt) {
01229     const unsigned int nEvt = evthdr.nevent; // evt->GetEntries();
01230 
01231     for (unsigned int i = 0; i < nEvt; i++){
01232       NtpSREvent* evtEntry = dynamic_cast<NtpSREvent*>( evt->At(i) );
01233       if ( !evtEntry ) continue;
01234 
01235       NtpTHEvent*  thevtEntry  = dynamic_cast<NtpTHEvent*>(thevt->At(i));
01236       NtpMCTruth*  mcEntry     = (thevtEntry)?
01237         dynamic_cast<NtpMCTruth*>( mc->At( (thevtEntry->neumc) ) ) : 0;
01238       if ( !mcEntry ) continue;
01239 
01240       const double eReco = evtEntry->ph.gev;
01241       const double eTrue = mcEntry->p4neu[3];
01242 
01243       if (eReco < eRecoMax || eTrue < eTrueMax ) continue;
01244 
01245       VHS::evtType evtt = VHS::GetEvtType(mcEntry->inu, mcEntry->iaction);
01246 
01247       if ( evtt == vhsUnknown ) continue;
01248 
01249       if (evtEntry->plane.end - evtEntry->plane.beg > cutPlanes) continue;
01250 
01251       if ( ( evtt == vhsNC    && NCevts    >= maxTrain ) ||
01252            ( evtt == vhsCCe   && eCCevts   >= maxTrain ) ||
01253            ( evtt == vhsCCmu  && muCCevts  >= maxTrain ) ||
01254            ( evtt == vhsCCtau && tauCCevts >= maxTrain )  )
01255         continue;
01256 
01257       // Gather our image and index sparse vectors
01258       std::vector<double> image;
01259       std::vector<int>    index;
01260       std::vector<double> theta;
01261       VHS::GetImage(evtEntry->stp, evtEntry->nstrip,
01262                     stp, nPlanes, nStrips, image, index, theta);
01263 
01264       // Break up into cases by event type
01265       switch (evtt) {
01266         case vhsNC :
01267           for (unsigned int ind = 0; ind < image.size(); ind++) {
01268             int    ii = index[ind];
01269             double xi = image[ind];
01270             int    ni = numNC[ii];
01271             
01272             avgNC[ii] = (avgNC[ii]*ni + xi   )/(ni+1);
01273             varNC[ii] = (varNC[ii]*ni + xi*xi)/(ni+1);
01274             numNC[ii]++;
01275           } // end outer for loop 
01276           NCevts++;
01277 
01278           break;
01279         case vhsCCe :
01280           for (unsigned int ind = 0; ind < image.size(); ind++) {
01281             int    ii = index[ind];
01282             double xi = image[ind];
01283             int    ni = numCCe[ii];
01284             
01285             avgCCe[ii] = (avgCCe[ii]*ni + xi   )/(ni+1);
01286             varCCe[ii] = (varCCe[ii]*ni + xi*xi)/(ni + 1);
01287             numCCe[ii]++;
01288           } // end outer for loop 
01289           eCCevts++;
01290 
01291           break;
01292         case vhsCCmu :
01293           for (unsigned int ind = 0; ind < image.size(); ind++) {
01294             int    ii = index[ind];
01295             double xi = image[ind];
01296             int    ni = numCCmu[ii];
01297             
01298             avgCCmu[ii] = (avgCCmu[ii]*ni + xi   )/(ni+1);
01299             varCCmu[ii] = (varCCmu[ii]*ni + xi*xi)/(ni+1);
01300             numCCmu[ii]++;
01301           } // end outer for loop 
01302           muCCevts++;
01303 
01304           break;
01305         case vhsCCtau :
01306           for (unsigned int ind = 0; ind < image.size(); ind++) {
01307             int    ii = index[ind];
01308             double xi = image[ind];
01309             int    ni = numCCtau[ii];
01310             
01311             avgCCtau[ii] = (avgCCtau[ii]*ni + xi   )/(ni+1);
01312             varCCtau[ii] = (varCCtau[ii]*ni + xi*xi)/(ni+1);
01313             numCCtau[ii]++;
01314           } // end outer for loop 
01315           tauCCevts++;
01316 
01317           break;
01318         case vhsUnknown :
01319           continue;
01320       } // end switch (evtType)
01321 
01322       nEvtUsed++;
01323 
01324     } // end for loop over events
01325   } // end if the ntpst arrays are valid
01326 
01327   return nEvtUsed;
01328 
01329 }

void VHS::TrainPost ( TFile *  outFile,
std::vector< double > &  avgNC,
std::vector< double > &  avgCCe,
std::vector< double > &  avgCCmu,
std::vector< double > &  avgCCtau,
std::vector< double > &  varNC,
std::vector< double > &  varCCe,
std::vector< double > &  varCCmu,
std::vector< double > &  varCCtau,
std::vector< int > &  numNC,
std::vector< int > &  numCCe,
std::vector< int > &  numCCmu,
std::vector< int > &  numCCtau,
int  NCevts,
int  eCCevts,
int  muCCevts,
int  tauCCevts,
const int  nPlanes = 20,
const int  nStrips = 20 
)

Definition at line 1333 of file VHS.cxx.

References WriteFile().

01352 {
01353   //
01354   // Post-processing of the variances, assuming input "variances"
01355   // are entries simply of <x^2>
01356   //
01357   // Translate averages & completed variances into TH2D and write
01358   // to file, outFile
01359   //
01360   // Call this only once at the end of the total training job
01361   //
01362 
01363   const unsigned int kImgSz   = avgNC.size();
01364 
01365   // Calculate the simple raw hit probabilities
01366   std::vector<double> pNChit   (kImgSz, 0.0);
01367   std::vector<double> pCCehit  (kImgSz, 0.0);
01368   std::vector<double> pCCmuhit (kImgSz, 0.0);
01369   std::vector<double> pCCtauhit(kImgSz, 0.0);
01370 
01371   for (unsigned int pix = 0; pix < kImgSz; pix++) {
01372     pNChit[pix]    = (NCevts > 0)?
01373       (double)numNC[pix]    / (double)NCevts    : 0.;
01374     pCCehit[pix]   = (eCCevts > 0)?
01375       (double)numCCe[pix]   / (double)eCCevts   : 0.;
01376     pCCmuhit[pix]  = (muCCevts > 0)?
01377       (double)numCCmu[pix]  / (double)muCCevts  : 0.;
01378     pCCtauhit[pix] = (tauCCevts > 0)?
01379       (double)numCCtau[pix] / (double)tauCCevts : 0.;
01380   }
01381 
01382   // Complete the variance calculation: var = <x^2> - <x>^2
01383   for (unsigned int ii = 0; ii < kImgSz; ii++) {
01384     varNC[ii]    -= avgNC[ii]    * avgNC[ii];
01385     varCCe[ii]   -= avgCCe[ii]   * avgCCe[ii];
01386     varCCmu[ii]  -= avgCCmu[ii]  * avgCCmu[ii];
01387     varCCtau[ii] -= avgCCtau[ii] * avgCCtau[ii];
01388 
01389   } // end strip loop
01390 
01391   VHS::WriteFile( outFile,
01392                   avgNC, avgCCe, avgCCmu, avgCCtau,
01393                   varNC, varCCe, varCCmu, varCCtau,
01394                   pNChit, pCCehit, pCCmuhit, pCCtauhit,
01395                   nPlanes, nStrips);
01396 
01397   return;
01398 
01399 }

void VHS::UnitVector ( std::vector< double > &  vec  ) 

Definition at line 1403 of file VHS.cxx.

References len.

Referenced by FillDiscriminants(), and FindMedian().

01404 {
01405   double len = 0.;
01406 
01407   for (unsigned int i = 0; i < vec.size(); i++)
01408     len += TMath::Power( vec[i], 2. );
01409 
01410   len = TMath::Sqrt( len );
01411 
01412   if (len <= 0.) return;
01413 
01414   for (unsigned int i = 0; i < vec.size(); i++)
01415     vec[i] = vec[i] / len ;
01416 
01417   return;
01418 
01419 }

int VHS::VecIndex ( int  ipln,
int  istp,
int  nPlanes 
)

Definition at line 1423 of file VHS.cxx.

Referenced by GetImage(), and ToVector().

01424 {
01425   // A particular choice of convention
01426 
01427   return (ipln + nPlanes * istp);
01428 
01429 }

void VHS::WriteFile ( TFile *  outFile,
std::vector< double > &  avgNC,
std::vector< double > &  avgCCe,
std::vector< double > &  avgCCmu,
std::vector< double > &  avgCCtau,
std::vector< double > &  varNC,
std::vector< double > &  varCCe,
std::vector< double > &  varCCmu,
std::vector< double > &  varCCtau,
std::vector< double > &  pNChit,
std::vector< double > &  pCCehit,
std::vector< double > &  pCCmuhit,
std::vector< double > &  pCCtauhit,
const int  nPlanes,
const int  nStrips 
)

Definition at line 1433 of file VHS.cxx.

References ToTH2D().

Referenced by TrainPost().

01448 {
01449   //
01450   // outFile must be open and valid
01451   //
01452   // Translate averages & completed variances into TH2D and write
01453   // to file, outFile
01454   //
01455   // Call this only once at the end of the total training job
01456   //
01457 
01458   // Prepare to write our training info to output file
01459   outFile->cd();
01460 
01461   const std::string imgXTitle = "Plane";
01462   const std::string imgYTitle = "Strip";
01463 
01464   // NC hit prob hist
01465   const std::string pNCTitle = "NC Hit Probabilities";
01466   const std::string pNCName  = "pNCHit";
01467   TH2D* UpNCHist; TH2D* VpNCHist;
01468   VHS::ToTH2D( pNChit,
01469                nPlanes, nStrips,
01470                pNCName.c_str(),
01471                pNCTitle.c_str(),
01472                UpNCHist, VpNCHist);
01473   UpNCHist->SetXTitle( imgXTitle.c_str() );
01474   UpNCHist->SetYTitle( imgYTitle.c_str() );
01475   UpNCHist->Write();
01476   VpNCHist->SetXTitle( imgXTitle.c_str() );
01477   VpNCHist->SetYTitle( imgYTitle.c_str() );
01478   VpNCHist->Write();
01479 
01480   // CCe hit prob hist
01481   const std::string pCCeTitle = "CCe Hit Probabilities";
01482   const std::string pCCeName  = "pCCeHit";
01483   TH2D* UpCCeHist; TH2D* VpCCeHist;
01484   VHS::ToTH2D( pCCehit,
01485                nPlanes, nStrips,
01486                pCCeName.c_str(),
01487                pCCeTitle.c_str(),
01488                UpCCeHist, VpCCeHist);
01489   UpCCeHist->SetXTitle( imgXTitle.c_str() );
01490   UpCCeHist->SetYTitle( imgYTitle.c_str() );
01491   UpCCeHist->Write();
01492   VpCCeHist->SetXTitle( imgXTitle.c_str() );
01493   VpCCeHist->SetYTitle( imgYTitle.c_str() );
01494   VpCCeHist->Write();
01495 
01496   // CCmu hit prob hist
01497   const std::string pCCmuTitle = "CC#mu Hit Probabilities";
01498   const std::string pCCmuName  = "pCCmuHit";
01499   TH2D* UpCCmuHist; TH2D* VpCCmuHist;
01500   VHS::ToTH2D( pCCmuhit,
01501                nPlanes, nStrips,
01502                pCCmuName.c_str(),
01503                pCCmuTitle.c_str(),
01504                UpCCmuHist, VpCCmuHist);
01505   UpCCmuHist->SetXTitle( imgXTitle.c_str() );
01506   UpCCmuHist->SetYTitle( imgYTitle.c_str() );
01507   UpCCmuHist->Write();
01508   VpCCmuHist->SetXTitle( imgXTitle.c_str() );
01509   VpCCmuHist->SetYTitle( imgYTitle.c_str() );
01510   VpCCmuHist->Write();
01511 
01512   // CCtau hit prob hist
01513   const std::string pCCtauTitle = "CC#tau Hit Probabilities";
01514   const std::string pCCtauName  = "pCCtauHit";
01515   TH2D* UpCCtauHist; TH2D* VpCCtauHist;
01516   VHS::ToTH2D( pCCtauhit,
01517                nPlanes, nStrips,
01518                pCCtauName.c_str(),
01519                pCCtauTitle.c_str(),
01520                UpCCtauHist, VpCCtauHist);
01521   UpCCtauHist->SetXTitle( imgXTitle.c_str() );
01522   UpCCtauHist->SetYTitle( imgYTitle.c_str() );
01523   UpCCtauHist->Write();
01524   VpCCtauHist->SetXTitle( imgXTitle.c_str() );
01525   VpCCtauHist->SetYTitle( imgYTitle.c_str() );
01526   VpCCtauHist->Write();
01527 
01528   // NC Avg hist
01529   const std::string avgNCTitle = "Average NC Image";
01530   const std::string avgNCName  = "AvgNCImg";
01531   TH2D* UavgNCHist; TH2D* VavgNCHist;
01532   VHS::ToTH2D( avgNC,
01533                nPlanes, nStrips,
01534                avgNCName.c_str(),
01535                avgNCTitle.c_str(),
01536                UavgNCHist, VavgNCHist);
01537   UavgNCHist->SetXTitle( imgXTitle.c_str() );
01538   UavgNCHist->SetYTitle( imgYTitle.c_str() );
01539   UavgNCHist->Write();
01540   VavgNCHist->SetXTitle( imgXTitle.c_str() );
01541   VavgNCHist->SetYTitle( imgYTitle.c_str() );
01542   VavgNCHist->Write();
01543 
01544   // CCe Avg hist
01545   const std::string avgCCeTitle = "Average CCe Image";
01546   const std::string avgCCeName  = "AvgCCeImg";
01547   TH2D* UavgCCeHist; TH2D* VavgCCeHist;
01548   VHS::ToTH2D( avgCCe,
01549                nPlanes, nStrips,
01550                avgCCeName.c_str(),
01551                avgCCeTitle.c_str(),
01552                UavgCCeHist, VavgCCeHist);
01553   UavgCCeHist->SetXTitle( imgXTitle.c_str() );
01554   UavgCCeHist->SetYTitle( imgYTitle.c_str() );
01555   UavgCCeHist->Write();
01556   VavgCCeHist->SetXTitle( imgXTitle.c_str() );
01557   VavgCCeHist->SetYTitle( imgYTitle.c_str() );
01558   VavgCCeHist->Write();
01559 
01560   // CCmu Avg hist
01561   const std::string avgCCmuTitle = "Average CC#mu Image";
01562   const std::string avgCCmuName  = "AvgCCmuImg";
01563   TH2D* UavgCCmuHist; TH2D* VavgCCmuHist;
01564   VHS::ToTH2D( avgCCmu,
01565                nPlanes, nStrips,
01566                avgCCmuName.c_str(),
01567                avgCCmuTitle.c_str(),
01568                UavgCCmuHist, VavgCCmuHist);
01569   UavgCCmuHist->SetXTitle( imgXTitle.c_str() );
01570   UavgCCmuHist->SetYTitle( imgYTitle.c_str() );
01571   UavgCCmuHist->Write();
01572   VavgCCmuHist->SetXTitle( imgXTitle.c_str() );
01573   VavgCCmuHist->SetYTitle( imgYTitle.c_str() );
01574   VavgCCmuHist->Write();
01575 
01576   // CCtau Avg hist
01577   const std::string avgCCtauTitle = "Average CC#tau Image";
01578   const std::string avgCCtauName  = "AvgCCtauImg";
01579   TH2D* UavgCCtauHist; TH2D* VavgCCtauHist;
01580   VHS::ToTH2D( avgCCtau,
01581                nPlanes, nStrips,
01582                avgCCtauName.c_str(),
01583                avgCCtauTitle.c_str(),
01584                UavgCCtauHist, VavgCCtauHist);
01585   UavgCCtauHist->SetXTitle( imgXTitle.c_str() );
01586   UavgCCtauHist->SetYTitle( imgYTitle.c_str() );
01587   UavgCCtauHist->Write();
01588   VavgCCtauHist->SetXTitle( imgXTitle.c_str() );
01589   VavgCCtauHist->SetYTitle( imgYTitle.c_str() );
01590   VavgCCtauHist->Write();
01591 
01592   // NC variance hist
01593   const std::string varNCTitle = "NC Variance";
01594   const std::string varNCName  = "NCVar";
01595   TH2D* UvarNCHist; TH2D* VvarNCHist;
01596   VHS::ToTH2D( varNC,
01597                nPlanes, nStrips,
01598                varNCName.c_str(),
01599                varNCTitle.c_str(),
01600                UvarNCHist, VvarNCHist);
01601   UvarNCHist->SetXTitle( imgXTitle.c_str() );
01602   UvarNCHist->SetYTitle( imgYTitle.c_str() );
01603   UvarNCHist->Write();
01604   VvarNCHist->SetXTitle( imgXTitle.c_str() );
01605   VvarNCHist->SetYTitle( imgYTitle.c_str() );
01606   VvarNCHist->Write();
01607 
01608   // CCe variance hist
01609   const std::string varCCeTitle = "CCe Variance";
01610   const std::string varCCeName  = "CCeVar";
01611   TH2D* UvarCCeHist; TH2D* VvarCCeHist;
01612   VHS::ToTH2D( varCCe,
01613                nPlanes, nStrips,
01614                varCCeName.c_str(),
01615                varCCeTitle.c_str(),
01616                UvarCCeHist, VvarCCeHist);
01617   UvarCCeHist->SetXTitle( imgXTitle.c_str() );
01618   UvarCCeHist->SetYTitle( imgYTitle.c_str() );
01619   UvarCCeHist->Write();
01620   VvarCCeHist->SetXTitle( imgXTitle.c_str() );
01621   VvarCCeHist->SetYTitle( imgYTitle.c_str() );
01622   VvarCCeHist->Write();
01623 
01624   // CCmu variance hist
01625   const std::string varCCmuTitle = "CCmu Variance";
01626   const std::string varCCmuName  = "CCmuVar";
01627   TH2D* UvarCCmuHist; TH2D* VvarCCmuHist;
01628   VHS::ToTH2D( varCCmu,
01629                nPlanes, nStrips,
01630                varCCmuName.c_str(),
01631                varCCmuTitle.c_str(),
01632                UvarCCmuHist, VvarCCmuHist);
01633   UvarCCmuHist->SetXTitle( imgXTitle.c_str() );
01634   UvarCCmuHist->SetYTitle( imgYTitle.c_str() );
01635   UvarCCmuHist->Write();
01636   VvarCCmuHist->SetXTitle( imgXTitle.c_str() );
01637   VvarCCmuHist->SetYTitle( imgYTitle.c_str() );
01638   VvarCCmuHist->Write();
01639 
01640   // CCtau variance hist
01641   const std::string varCCtauTitle = "CCtau Variance";
01642   const std::string varCCtauName  = "CCtauVar";
01643   TH2D* UvarCCtauHist; TH2D* VvarCCtauHist;
01644   VHS::ToTH2D( varCCtau,
01645                nPlanes, nStrips,
01646                varCCtauName.c_str(),
01647                varCCtauTitle.c_str(),
01648                UvarCCtauHist, VvarCCtauHist);
01649   UvarCCtauHist->SetXTitle( imgXTitle.c_str() );
01650   UvarCCtauHist->SetYTitle( imgYTitle.c_str() );
01651   UvarCCtauHist->Write();
01652   VvarCCtauHist->SetXTitle( imgXTitle.c_str() );
01653   VvarCCtauHist->SetYTitle( imgYTitle.c_str() );
01654   VvarCCtauHist->Write();
01655 
01656   return;
01657 
01658 }


Generated on 3 Dec 2018 for loon by  doxygen 1.6.1