ParticleFinder Class Reference

#include <ParticleFinder.h>

Inheritance diagram for ParticleFinder:
JobCModule

List of all members.

Public Member Functions

 ParticleFinder ()
 ~ParticleFinder ()
void BeginJob ()
void EndJob ()
JobCResult Reco (MomNavigator *mom)
const RegistryDefaultConfig () const
void Config (const Registry &r)
void Reset ()
void FillPot (ParticleBeamMon *bmon)

Private Member Functions

double OscillationProb (TF1 *f2, int ntype, double NuE, double sinth23, double sin2th13)
void CheckContainment (ParticleObjectHolder *p, double front_z[2][nfront], double front_t[2][nfront])
int IsInsideNearFiducial_Nue_Standard (double x, double y, double z, bool)
int IsInsideFarFiducial_Nue_Standard (double x, double y, double z, bool isMC)

Private Attributes

Finder finder
std::string kPOTTreeName
int lastrun
ParticleBeamMonAna bma
int recoCnt
std::string outObjectName
int DoMRCC
double pecutlevel
double MIPperPE
double MEUperGeV

Static Private Attributes

static TF1 * osceq = 0
static SKZPWeightCalculatorskzpCalc = 0

Detailed Description

Definition at line 31 of file ParticleFinder.h.


Constructor & Destructor Documentation

ParticleFinder::ParticleFinder (  ) 

(Document me!)

Definition at line 70 of file ParticleFinder.cxx.

References DoMRCC, lastrun, MIPperPE, osceq, pecutlevel, recoCnt, and skzpCalc.

00071 {
00072         lastrun=-1;
00073 
00074         recoCnt=0;
00075         pecutlevel=0;
00076         MIPperPE=60.;
00077         DoMRCC=0;
00078 
00079 //      if(!pot)pot = new POT();
00080         if(!skzpCalc)skzpCalc=new SKZPWeightCalculator("DetXs",true);
00081   
00082   
00083 
00084         if(!osceq)osceq = new TF1("f2","sin(1.267*[0]*[1]/x)*sin(1.267*[0]*[1]/x)",0.,120.);
00085 
00086  
00087 
00088 
00089         
00093 }

ParticleFinder::~ParticleFinder (  ) 

(Document me!)

Definition at line 96 of file ParticleFinder.cxx.

References osceq, and skzpCalc.

00097 {
00101 
00102   
00103   //if(pot)delete pot;
00104   //pot=0;
00105   if(skzpCalc)delete skzpCalc;
00106   skzpCalc=0;
00107   if(osceq)delete osceq;osceq=0;
00108   //kb =0;
00109 }


Member Function Documentation

void ParticleFinder::BeginJob ( void   )  [virtual]

Implement for notification of begin of job

(Document me!)

Reimplemented from JobCModule.

Definition at line 113 of file ParticleFinder.cxx.

00114 {
00118  
00119   
00120 }

void ParticleFinder::CheckContainment ( ParticleObjectHolder p,
double  front_z[2][nfront],
double  front_t[2][nfront] 
) [private]

Definition at line 690 of file ParticleFinder.cxx.

References ParticleEvent::contained, ParticleObjectHolder::event, VldContext::GetDetector(), RecRecordImp< T >::GetHeader(), RecHeader::GetVldContext(), Detector::kNear, and tc.

Referenced by Reco().

00691 {
00692 
00693         if(p->GetHeader().GetVldContext().GetDetector()!=Detector::kNear)       
00694         {
00695                 p->event.contained=1;
00696                 return;
00697         }       
00698         
00699         //we are in the near... make sure that the event starts in the fully instrumented region........
00700 
00701         //here we use a fidual volume like check....
00702         //actual dimensions taken from viewing strip.tpos vs strip.z is as folows:
00703         // u view:  -0.2<t<2.3  0<z<7
00704         // v view:  -2.3<t<0.2  0<z<7
00705 
00706 
00707         double SuperModule1Beg = 0.1;//1.01080;  //Data and MC values (according to DataUtil/infid.h on 10/02/07
00708         double SuperModule1End = 7;//4.99059;
00709         double radialInner = 0;
00710         double radialOuter = 1.;//0.8;
00711         
00712         double uCenter = 1.1513;
00713         double vCenter = -0.9537;
00714         
00715         
00716         int contained=1;
00717         for(int view=0;view<2;view++)
00718         for(int c=0;c<nfront;c++)
00719         {
00720         
00721                 if(front_z[view][c]>=10000)continue;//not used.... its a short event
00722                 
00723                 
00724                 if(front_z[view][c]<SuperModule1Beg || front_z[view][c]>SuperModule1End)
00725                 {
00726                         contained=0;
00727                         break;
00728                 }
00729                 
00730                 double tc=0;
00731                 
00732                 if(view==0)
00733                 {
00734                         tc=front_t[view][c]-uCenter;
00735                 }else{
00736                         tc=front_t[view][c]-vCenter;
00737                 }
00738                 
00739                 tc=fabs(tc);//because its a radial cut and 0 is the center...
00740         
00741                 if(tc < radialInner || tc > radialOuter)
00742                 {
00743                         //printf("failing contained view %d st %f t %f z %f\n",view, tc,front_t[view][c], front_z[view][c]);
00744                         contained=0;
00745                         break;
00746                 }
00747                 
00748         }
00749 
00750         p->event.contained=contained;
00751 
00753 /*
00754 Double_t SuperModule1Beg = 1.01080;  //Data and MC values (according to DataUtil/infid.h on 10/02/07
00755   Double_t SuperModule1End = 4.99059;
00756 
00757     fX = 1./sqrt(2.0)*(fU - fV);
00758     fY = 1./sqrt(2.0)*(fU + fV);
00759 
00760         fu = 1./sqrt(2.)*(fx+fy);
00761         fv = 1./sqrt(2.)*(fy-fx);
00762 
00763   Double_t radialInner = 0;
00764   Double_t radialOuter = 0.8;
00765   Double_t xCenter = 1.4885;
00766   Double_t yCenter = 0.1397;
00767 
00768         uCenter = 1.1513;
00769         vCenter = -0.9537;
00770 
00771 
00772 */
00774         
00775                 
00776 
00777 }

void ParticleFinder::Config ( const Registry r  )  [virtual]

Return the actual configuration. If your module directly pulls its configuration from the fConfig Registry, you don't need to override this. Override if you have local config variables.

Configure the module given the Registry r

Reimplemented from JobCModule.

Definition at line 821 of file ParticleFinder.cxx.

References DoMRCC, Registry::Get(), MEUperGeV, MIPperPE, outObjectName, and pecutlevel.

00822 {
00826   
00827       const char* tmps;
00828         int tmpi;
00829         double tmpd;
00830   //  if(r.Get("POTTreeFileName",tmps)){kPOTTreeName=tmps;}
00831         if(r.Get("OutObjectName",tmps))
00832         {
00833                 outObjectName=tmps;
00834         //      cout << "Setting Particle Finder output object name to "<<outObjectName<<"\n";
00835         }
00836         if(r.Get("MEUperGeV",tmpd)){MEUperGeV=tmpd;}
00837         if(r.Get("DoMRCC",tmpi)){DoMRCC=tmpi;}
00838         if(r.Get("PECutLevel",tmpd)){pecutlevel=tmpd;}
00839         if(r.Get("MIPperPE",tmpd)){MIPperPE=tmpd;}
00840 
00841 
00842 }

const Registry & ParticleFinder::DefaultConfig ( void   )  const [virtual]

Get the default configuration registry. This should normally be overridden. One useful idiom is to implement it like:

const Registry& MyModule::DefaultConfig() const { static Registry cfg; // never is destroyed if (cfg.Size()) return cfg; // already filled it // set defaults: cfg.Set("TheAnswer",42); cfg.Set("Units","unknown"); return cfg; }

Supply the default configuration for the module

Reimplemented from JobCModule.

Definition at line 785 of file ParticleFinder.cxx.

References JobCModule::GetName(), Registry::LockValues(), Registry::Set(), and Registry::UnLockValues().

00786 {
00790   static Registry r; // Default configuration for module
00791 
00792   // Set name of config
00793   std::string name = this->GetName();
00794   name += ".config.default";
00795   r.SetName(name.c_str());
00796 
00797   //r.Set("POTTreeFileName","pottree.root");
00798   
00799 
00800   // Set values in configuration
00801   r.UnLockValues();
00802 
00803   r.Set("DoMRCC",0);
00804   r.Set("OutObjectName","Normal");
00805   r.Set("PECutLevel",2);
00806   r.Set("MIPperPE",0.1);
00807   r.Set("MEUperGeV",25.);
00808   r.LockValues();
00809 
00810 
00811   
00812 
00813 
00814 
00815 
00816   return r;
00817 }

void ParticleFinder::EndJob (  )  [virtual]

Implement for notification of end of job

(Document me!)

Reimplemented from JobCModule.

Definition at line 124 of file ParticleFinder.cxx.

00125 {
00129 
00130 /*
00131   //MSG("ParticleFinder",Msg::kInfo)
00132   cout<<"Number of POT in this job: "<<pot->pot<<endl;
00133         
00134    TDirectory *savedir = gDirectory;
00135         
00136    TFile *fpf = dynamic_cast<TFile *>(gROOT->GetListOfFiles()->FindObject(kPOTTreeName.c_str()));
00137    if(fpf){
00138      fpf->cd();
00139      TTree *pottree = new TTree("pottree","pottree");
00140      pottree->Branch("ParticlePOT",&pot);
00141      pottree->Fill();
00142      pottree->Write();
00143      savedir->cd();
00144    }
00145    else{
00146      MSG("ParticleFinder",Msg::kError)<<"Could not find the file to write the pottree to, there will be no pottree"<<endl;
00147    }
00148 
00149 */
00150 
00151 }

void ParticleFinder::FillPot ( ParticleBeamMon bmon  ) 

Definition at line 856 of file ParticleFinder.cxx.

References RecRecordImp< T >::GetHeader(), RecDataHeader::GetRun(), VldContext::GetSimFlag(), RecHeader::GetVldContext(), SimFlag::kMC, and lastrun.

Referenced by Reco().

00857 {
00858                 int ismc=bmon->GetHeader().GetVldContext().GetSimFlag() == SimFlag::kMC;
00859                 VldContext evt_vldc = bmon->GetHeader().GetVldContext();
00860 
00861 
00862                 if(!ismc)
00863                 {
00864                         //pot->pot+=bmon->GetPot();
00865                 }else{
00866                 
00867                         std::string relName = bmon->GetTitle();
00868 
00869                         std::string mcinfo = "";
00870                         if(bmon->GetHeader().GetVldContext().GetSimFlag() == SimFlag::kMC){
00871                         mcinfo = "Daikon";
00872         //              string temp = mchdr.geninfo.codename;
00873         //              if(temp.size() != 0){   mcinfo = temp;  }
00874                         }
00875                                         
00876                         //ReleaseType::Release_t rel = ReleaseType::MakeReleaseType(relName, mcinfo);
00877                 
00878                         
00879                         //near for every snarl
00880 /*                      if(evt_vldc.GetDetector()==Detector::kNear){
00881               double thispot = MCInfo::GetMCPoT(Detector::kNear, bmon->beamtype, rel);
00882               pot->pot+=thispot;
00883            // printf("near pot %f bt %d r %d\n",thispot,bmon->beamtype,rel);
00884             }
00885 */                      
00886                         //far only once per run
00887 /*                      if(evt_vldc.GetDetector()==Detector::kFar && 
00888                 bmon->GetHeader().GetRun()!=lastrun){
00889                         double thispot=6.5e8;//
00890                         //dogwood0 isn't showing up as daikon....
00891                         //MCInfo::GetMCPoT(Detector::kFar, bmon->beamtype, rel);
00892                         pot->pot+=thispot;
00893                         //printf("far pot %f bt %d r %d\n",thispot,bmon->beamtype,rel);
00894                 }
00895   */  
00896                 }
00897                 
00898                 
00899                 
00900 /*              pot->nsnarls++;
00901                 if(bmon->GetHeader().GetRun()!=lastrun)
00902                 {
00903                         pot->nruns++;
00904                         pot->beamtype = bmon->beamtype;
00905                 }
00906 */              lastrun=bmon->GetHeader().GetRun();
00907                 
00908 
00909 }

int ParticleFinder::IsInsideFarFiducial_Nue_Standard ( double  x,
double  y,
double  z,
bool  isMC 
) [private]

Definition at line 995 of file ParticleFinder.cxx.

Referenced by Reco().

00995                                                                                            {
00996 
00997 //cout << "vtx "<< x <<" "<<y<<" "<<z<<" "<<isMC<<endl;
00998 
00999   Double_t SuperModule1Beg =  0.49080;   // These are data values
01000   Double_t SuperModule2Beg = 16.27110;
01001   Double_t SuperModule1End = 14.29300;
01002   Double_t SuperModule2End = 27.98270;
01003 
01004   if(isMC){
01005     SuperModule1Beg =  0.47692;   // These are mc values
01006     SuperModule2Beg = 16.26470;
01007     SuperModule1End = 14.27860;
01008     SuperModule2End = 27.97240;
01009   }
01010 
01011   Double_t radialInner = 0.50;
01012   Double_t radialOuter = TMath::Sqrt(14.0);
01013   Bool_t zContained = false;
01014   Bool_t xyContained = false;
01015 
01016   Double_t r = TMath::Sqrt(x*x + y*y);
01017 
01018   if( (z >= SuperModule1Beg && z <=SuperModule1End) ||
01019       (z >= SuperModule2Beg && z <=SuperModule2End) )
01020      zContained = true;
01021 
01022   if( r >= radialInner && r <= radialOuter)
01023      xyContained = true;
01024 
01025   Int_t retVal = 0;
01026   if(zContained && xyContained) retVal = 1;
01027   if(!zContained) retVal = -1;
01028   if(!xyContained) retVal -= 2;
01029 
01030   return retVal;  //  1 contained, -1 out of bounds z
01031                   //  -2 oob xy, -3 oob both
01032 }

int ParticleFinder::IsInsideNearFiducial_Nue_Standard ( double  x,
double  y,
double  z,
bool   
) [private]

Definition at line 967 of file ParticleFinder.cxx.

Referenced by Reco().

00968 {
00969   Double_t SuperModule1Beg = 1.01080;  //Data and MC values (according to DataUtil/infid.h on 10/02/07
00970   Double_t SuperModule1End = 4.99059;
00971 
00972   Double_t radialInner = 0;
00973   Double_t radialOuter = 0.8;
00974   Double_t xCenter = 1.4885;
00975   Double_t yCenter = 0.1397;
00976 
00977   Bool_t zContained = false;
00978   Bool_t xyContained = false;
00979 
00980   Double_t r = TMath::Sqrt((x-xCenter)*(x-xCenter) + (y-yCenter)*(y-yCenter));
00981 
00982   if( z >= SuperModule1Beg && z <=SuperModule1End)
00983      zContained = true;
00984   if( r >= radialInner && r <= radialOuter)
00985      xyContained = true;
00986 
00987   Int_t retVal = 0;
00988   if(zContained && xyContained) retVal = 1;
00989   if(!zContained) retVal = -1;
00990   if(!xyContained) retVal -= 2;
00991 
00992   return retVal;
00993 }

double ParticleFinder::OscillationProb ( TF1 *  f2,
int  ntype,
double  NuE,
double  sinth23,
double  sin2th13 
) [private]

Definition at line 919 of file ParticleFinder.cxx.

References UtilMisc::OscProb().

Referenced by Reco().

00919                                                                                                       {
00920 
00921  // sinth23 : sine square theta 23
00922  // sin2th13 : sine square 2 theta 13
00923  //
00924  double OscProb = 0 ;
00925  double NumuToNutau ;
00926  double NumuToNue ;
00927  double NueSurvival ;
00928  double NumuSurvival ;
00929  double NueToNutau  ;
00930  double NueToNumu;
00931 
00932  NumuToNutau = 4.*sinth23*(1.-sinth23)*pow(1-sin2th13/4,2) ;
00933  NumuToNutau *= f2->Eval(TMath::Abs(NuE)) ;
00934 
00935  NumuToNue = sinth23*sin2th13*f2->Eval(TMath::Abs(NuE)) ;
00936  
00937 //printf("ting p = %f  sinth23 %f sin2th13 %f * l/e %f\n",NumuToNue,sinth23,sin2th13,f2->Eval(TMath::Abs(NuE))) ;
00938 
00939 //printf("ting l/e term %f\n",f2->Eval(TMath::Abs(NuE))) ;
00940 
00941  NueSurvival = 1.- sin2th13*f2->Eval(TMath::Abs(NuE)) ;
00942  //NueSurvival = 1.;
00943 
00944  NumuSurvival = 1. - NumuToNutau - NumuToNue ;
00945  //NumuSurvival = 1.;
00946 
00947  NueToNutau = (1.-sinth23)*sin2th13*f2->Eval(TMath::Abs(NuE)) ;
00948 
00949  NueToNumu = NumuToNue;
00950 
00951  if (ntype==0) OscProb = 1 ;
00952  else if (ntype==1) OscProb = NumuSurvival ;
00953  else if (ntype==2) OscProb = NumuToNue ;
00954  else if (ntype==3) OscProb = NumuToNutau ;
00955  else if (ntype==4) OscProb = NueSurvival ;
00956  else if (ntype==5) OscProb = NueToNutau ;
00957  else if (ntype==6) OscProb = NueToNumu;
00958 //  cout<<"Period = "<<f2->Eval(TMath::Abs(NuE))<<endl ;
00959 //  cout<<"Oscillation probability = "<<OscProb<<endl ;
00960 
00961  return OscProb ;
00962 } 

JobCResult ParticleFinder::Reco ( MomNavigator mom  )  [virtual]

Implement this for read-write access to the MomNavigator

Reimplemented from JobCModule.

Definition at line 155 of file ParticleFinder.cxx.

References MuELoss::a, Finder::AddStrip(), MomNavigator::AdoptFragment(), ParticleBeamMonAna::ana(), ParticleBeamMon::beamtype, NtpMREvent::best_complete, NtpMREvent::best_event, bma, CheckContainment(), NtpMCStdHep::child, det, DoMRCC, exit(), FillPot(), finder, ParticleBeamMonAna::fname, SKZPWeightCalculator::GetBeamWeight(), VldContext::GetDetector(), RecDataHeader::GetErrorCode(), NtpStRecord::GetEvents(), MomNavigator::GetFragment(), MomNavigator::GetFragmentList(), RecRecordImp< T >::GetHeader(), NtpStRecord::GetMCStdHeps(), NtpStRecord::GetMCTruths(), SntpHelpers::GetMREvent(), NtpStRecord::GetRelease(), RecPhysicsHeader::GetRemoteSpillType(), DbiResultPtr< T >::GetRow(), RecDataHeader::GetRun(), RecDataHeader::GetRunType(), CalMIPCalibration::GetScale(), VldContext::GetSimFlag(), RecPhysicsHeader::GetSnarl(), Finder::GetStrips(), NtpStRecord::GetStrips(), RecDataHeader::GetSubRun(), NtpStRecord::GetTHEvents(), RecPhysicsHeader::GetTimeFrame(), RecPhysicsHeader::GetTrigSrc(), RecHeader::GetVldContext(), NtpMCStdHep::index, Managed::HitManager::InsertHit(), ReleaseType::IsBirch(), IsInsideFarFiducial_Nue_Standard(), IsInsideNearFiducial_Nue_Standard(), Msg::kDebug, JobCResult::kFailed, Detector::kFar, SimFlag::kMC, Detector::kNear, JobCResult::kPassed, kPOTTreeName, Msg::kWarning, MEUperGeV, MIPperPE, NtpMRRecord::mrhdr, MSG, n, DbiResultPtr< T >::NewQuery(), nfront, NtpMRSummary::nmrevt, NtpMREvent::orig_event, osceq, OscillationProb(), outObjectName, NtpMCStdHep::parent, pecutlevel, Finder::Process(), recoCnt, Finder::Reset(), Finder::SetClusterManagerU(), Finder::SetClusterManagerV(), Managed::ClusterManager::SetClusterSaver(), Managed::ClusterManager::SetHitManager(), Finder::SetMEUperGeV(), Finder::SetMRCC(), Finder::SetPOH(), Finder::SetTrueVertex(), skzpCalc, BeamType::ToZarko(), NtpMCFluxInfo::tptype, NtpMCFluxInfo::tpx, NtpMCFluxInfo::tpy, NtpMCFluxInfo::tpz, Finder::vtx_u, Finder::vtx_v, and Finder::vtx_z.

00156 {
00157 
00158 
00159    bool foundMR=false;
00160    bool foundST=false;
00161 
00162   NtpStRecord *str=dynamic_cast<NtpStRecord *>(mom->GetFragment("NtpStRecord","Primary"));
00163   NtpMRRecord *mr = 0;
00164   NtpStRecord *oldstr = 0;
00165 
00166   VldContext vc;  
00167    
00168   if(str){
00169         foundST=true;   
00170         vc=str->GetHeader().GetVldContext();
00171         ReleaseType::Release_t release = str->GetRelease();
00172         //check for MR:
00173         mr =  dynamic_cast<NtpMRRecord *>(mom->GetFragment("NtpMRRecord"));
00174         if(mr){
00175                 if(ReleaseType::IsBirch(release)) {
00176                                 printf("please don't run with birch (particle finder)\n");exit(1);
00177                                 //oldstr=dynamic_cast<NtpStRecord *>(mom->GetFragment("NtpStRecord",
00178                                 //                          "NtpStRecordOld")); 
00179                 } else {
00180                         oldstr=str;
00181                         str=dynamic_cast<NtpStRecord *>(mom->GetFragment("NtpStRecord",
00182                                                           "MuonRemoved"));
00183                 }
00184                 if(oldstr) foundMR = true;
00185         }
00186   }
00187 
00188   if(!foundST && !foundMR)
00189   {
00190       MSG("ParticleFinder",Msg::kWarning)<<"Got Nothing from mom"<<endl;
00191   
00192       return JobCResult::kFailed;
00193 
00194   }
00195 
00196   //std::vector<ParticleObjectHolder *>* hf;
00197   std::vector<TObject* > hft =( mom->GetFragmentList("ParticleObjectHolder"));
00198 
00199 MSG("ParticleFinder",Msg::kDebug) << "Snarl " << str->GetHeader().GetSnarl() << endl;
00200 
00201   int make_new_holder=1;
00202 
00203 
00204   if(hft.size()>0)
00205   {
00206       make_new_holder=0;
00207   }
00208  
00209 
00210 MSG("ParticleFinder",Msg::kDebug) << "size of found array " << hft.size()<<endl; 
00211 
00212   int ismc=str->GetHeader().GetVldContext().GetSimFlag() == SimFlag::kMC;
00213 
00214 
00215   std::vector<const NtpSREvent* > evts = str->GetEvents();
00216   std::vector<const NtpSRStrip* > stp = str->GetStrips();
00217 
00218         std::vector<const NtpMCStdHep* > stdhep;
00219     std::vector<const NtpTHEvent* > thevt;
00220         std::vector<const NtpMCTruth* > mc;
00221 
00222         if(ismc)
00223         {
00224                 stdhep = oldstr ? oldstr->GetMCStdHeps() : str->GetMCStdHeps();
00225                 thevt = oldstr ? oldstr->GetTHEvents() : str->GetTHEvents();
00226                 mc = oldstr ? oldstr->GetMCTruths() : str->GetMCTruths();
00227         }
00228         
00229         
00230   MSG("ParticleFinder",Msg::kDebug) << "Events " << evts.size() <<endl;
00231   if(evts.size()<1)
00232         {
00233 
00234 /*              if(make_new_holder)
00235                 {
00236                                 RecCandHeader ntphdr(str->GetHeader().GetVldContext(),str->GetHeader().GetRun(),
00237                                         str->GetHeader().GetSubRun(),str->GetHeader().GetRunType(),str->GetHeader().GetErrorCode(),
00238                                         str->GetHeader().GetSnarl(),str->GetHeader().GetTrigSrc(),str->GetHeader().GetTimeFrame(),
00239                                         str->GetHeader().GetRemoteSpillType(),-1);
00240                 
00241                 
00242                         ParticleObjectHolder * h = new ParticleObjectHolder(ntphdr); 
00243                         
00244                 
00245                         mom->AdoptFragment(h);
00246                 }
00247 */
00248                 return JobCResult::kPassed;  
00249         }
00250 
00251         ParticleBeamMon *bmon=0;
00252         MSG("ParticleFinder",Msg::kDebug) << "!!! "<<evts.size()<<" events, "<<  hft.size()<<" currently there"<<endl;
00253 
00254 for(unsigned int ievt=0;ievt<evts.size();ievt++)
00255 {
00256 
00257 
00258         int bestmatch=ievt;
00259 
00260         if(mr)
00261         {
00262 
00263 
00264            bestmatch=-1;
00265            //Int_t best_rmmu = i;
00266            Float_t best_com = 0;
00267            for(int xi=0;xi<mr->mrhdr.nmrevt;xi++){        
00268              NtpMREvent *ev = SntpHelpers::GetMREvent(xi,mr);
00269              if(ev && ev->best_event==(int)ievt && ev->best_complete>best_com) {
00270                 best_com = ev->best_complete;
00271               // best_rmmu = ev->orig_event;
00272                 bestmatch = ev->orig_event;       
00273         //             cout<<"found mr match "<<i<<" "<<ev->best_event<<" "<<xi<<endl;
00274              }
00275            }
00276         }
00277 
00278 
00279 
00280 //int ievt=0;
00281 
00282         
00283 //      if(ievt!=4 ||str->GetHeader().GetSnarl()!=2800034)continue;
00284 //if(ievt!=0)continue;
00285 
00286         ParticleObjectHolder *n=0;
00287 
00288 
00289     //    if(make_new_holder || hft.size() < ievt+1)
00290      //   {
00291                         RecCandHeader ntphdr(str->GetHeader().GetVldContext(),str->GetHeader().GetRun(),
00292                                         str->GetHeader().GetSubRun(),str->GetHeader().GetRunType(),str->GetHeader().GetErrorCode(),
00293                                         str->GetHeader().GetSnarl(),str->GetHeader().GetTrigSrc(),str->GetHeader().GetTimeFrame(),
00294                                         str->GetHeader().GetRemoteSpillType(),ievt);
00295                 n = new ParticleObjectHolder(ntphdr);
00296 
00297                 
00298     //    }else{
00299      //           n=dynamic_cast<ParticleObjectHolder *>(hft[ievt]);
00300     //    }
00301 
00302 
00303         MSG("ParticleFinder",Msg::kDebug)<<"Generating event for run "<<str->GetHeader().GetRun()<<" subrun "<<str->GetHeader().GetSubRun()<<" snarl "<<str->GetHeader().GetSnarl()<<" event "<<ievt<<" matchingevt "<<bestmatch<<"\n";
00304         int mc_idx=-1;
00305         if(ismc && bestmatch>-1)mc_idx=thevt[bestmatch]->neumc;
00306                 
00307 
00308 
00309 if(ievt==0 || !bmon)
00310 {
00311 
00312         bmon = new ParticleBeamMon(ntphdr);
00313                 mom->AdoptFragment(bmon);
00314         
00315 
00316                 bma.fname=kPOTTreeName;
00317                 bma.ana(n,bmon);
00318 
00319                 
00320                 FillPot(bmon);
00321 }
00322 
00323 MSG("ParticleFinder",Msg::kDebug) << "partices in current evt "<< ievt << " array " << n->particles->GetEntries() << endl; 
00324 
00325         if(n->particles->GetEntries()>0)return JobCResult::kPassed;//continue; //dont rerun it
00326 
00327 
00328         finder.Reset();
00329         finder.SetMRCC(DoMRCC);
00330         finder.SetPOH(n);
00331         finder.SetMEUperGeV(MEUperGeV);
00332 
00333 
00334 
00335 
00336 
00337 
00338 
00339         if(ismc)
00340         {
00341         int beststdhep = -1;
00342         if(bestmatch>-1)beststdhep = thevt[bestmatch]->neustdhep;
00343         
00344         if(beststdhep>-1)//we may not have a match!
00345         {       
00346                 double vx = stdhep[beststdhep]->vtx[0];
00347                 double vy = stdhep[beststdhep]->vtx[1];
00348                 double vz = stdhep[beststdhep]->vtx[2];
00349         
00350                 finder.SetTrueVertex((vx+vy)/sqrt(2.0),(vy-vx)/sqrt(2.0),vz);
00351 
00352                 n->mctrue.vtx_u=(vx+vy)/sqrt(2.0);
00353                 n->mctrue.vtx_v=(-vx+vy)/sqrt(2.0);
00354                 n->mctrue.vtx_z=vz;
00355         
00356                 //want to save stdhep array...
00357                 int thismc = stdhep[beststdhep]->mc;
00358 
00359                 for(unsigned int i=beststdhep;i<stdhep.size() && stdhep[i]->mc==thismc;i++)
00360                 {
00361                 //      new(&(n->mctrue.stdhep[i-beststdhep])) NtpMCStdHep();
00362                 //      *((NtpMCStdHep *)(n->mctrue.stdhep->AddrAt(i-beststdhep)))=*(stdhep[i]);
00363                         
00364                 //      NtpMCStdHep *me = new NtpMCStdHep();
00365                 //      n->mctrue.stdhep->AddAtAndExpand(me,i-beststdhep);
00366                 //      *me = *(stdhep[i]);
00367         
00368                         NtpMCStdHep a=*(stdhep[i]);
00369         
00370                         a.parent[0]=a.parent[0]-beststdhep;
00371                         a.parent[0]=a.parent[1]-beststdhep;
00372                         a.child[0]=a.child[0]-beststdhep;
00373                         a.child[0]=a.child[1]-beststdhep;
00374                         a.index=i;
00375                         n->mctrue.stdhep.push_back(a);
00376         
00377         
00378                         if(mc_idx>-1)n->mctrue.visible_energy=mc[mc_idx]->tphu+mc[mc_idx]->tphv;
00379                 }
00380         //      printf("saved %d stdheps\n",n->mctrue.stdhep.size());
00381         
00382         }
00383         }
00384 
00385 
00386 
00387 
00388         
00389         n->event.sr_vtx_u=evts[0]->vtx.u;
00390         n->event.sr_vtx_v=evts[0]->vtx.v;
00391         n->event.sr_vtx_z=evts[0]->vtx.z;       
00392         
00393         //set this to what we are actually filling based on sigcors
00394         n->event.visenergy=-1;//evts[ievt]->ph.mip;
00395                 
00396         if(ismc && mc_idx>-1)
00397         {
00398 
00399                 n->mctrue.iaction=mc[mc_idx]->iaction;
00400                 n->mctrue.inu=mc[mc_idx]->inu;
00401                 n->mctrue.iresonance=mc[mc_idx]->iresonance;
00402                 n->mctrue.nuenergy=mc[mc_idx]->p4neu[3];
00403 
00404                 n->mctrue.inunoosc=mc[mc_idx]->inunoosc;
00405 
00406                 n->mctrue.flux=mc[mc_idx]->flux;
00407                 
00408         /*      n->mctrue.osc_dm2=0.0024;
00409                 n->mctrue.osc_L=735.; 
00410                 n->mctrue.osc_sinth23=0.5;
00411                 n->mctrue.osc_sin2th13=0.15;
00412         */      
00413                 n->mctrue.type=(n->mctrue.iaction>0)*(abs(n->mctrue.inu)/2-5);
00414                 //and for beam ves
00415                 if(n->mctrue.inu==n->mctrue.inunoosc && TMath::Abs(n->mctrue.inu)==12)n->mctrue.type=4;
00416                 
00417                 
00418         }
00419         
00420         
00421         
00423         
00424         if(ismc)
00425         {
00426         VldContext evt_vldc = n->GetHeader().GetVldContext();
00427 
00428 
00429         int det=evt_vldc.GetDetector(); 
00430         BeamType::BeamType_t beam=bmon->beamtype;
00431     NtpMCFluxInfo fi= n->mctrue.flux;
00432         double pt = sqrt(fi.tpx*fi.tpx+fi.tpy*fi.tpy);
00433         double pz = 1.*fi.tpz;
00434         int tptype = fi.tptype;
00435         int zbeam = BeamType::ToZarko(beam);
00436         n->mctrue.totbeamweight = skzpCalc->GetBeamWeight(det,zbeam,tptype,pt,pz,n->mctrue.nuenergy,n->mctrue.inu);
00437         }else{
00438                 n->mctrue.totbeamweight = 1;
00439         }
00441 
00442 
00444         if(ismc && n->GetHeader().GetVldContext().GetDetector()==Detector::kFar)
00445         {
00446         int otype=0;
00447         if(n->mctrue.iaction==1)
00448         {
00449                 if(abs(n->mctrue.inu)==12)
00450                 {
00451                         if(abs(n->mctrue.inunoosc)==12) otype=4;
00452                         else otype=2;
00453                 }       
00454                 if(abs(n->mctrue.inu)==14)
00455                 {
00456                         if(abs(n->mctrue.inunoosc)==12) otype=6;
00457                         else otype=1;
00458                 }       
00459                 if(abs(n->mctrue.inu)==16)
00460                 {
00461                         if(abs(n->mctrue.inunoosc)==12) otype=5;
00462                         if(abs(n->mctrue.inunoosc)==14) otype=3;
00463                 }       
00464         }
00465 
00466 
00467         n->mctrue.osc_L=735;
00468         n->mctrue.osc_dm2=0.0024;
00469         n->mctrue.osc_sinth23=sin(3.141592/4.0);
00470         n->mctrue.osc_sin2th13=0.15;
00471         osceq->SetParameters(n->mctrue.osc_dm2,n->mctrue.osc_L);
00472         n->mctrue.oscprob=OscillationProb(osceq, otype, n->mctrue.nuenergy, n->mctrue.osc_sinth23* n->mctrue.osc_sinth23, n->mctrue.osc_sin2th13); 
00473 
00474 //      printf("type %d oscprob %f e %f\n",n->mctrue.type,n->mctrue.oscprob, n->mctrue.nuenergy);
00475 
00476         }
00478         
00479 
00480 
00481         Managed::HitManager hitmanager_u;
00482         Managed::ClusterManager clustermanager_u;
00483         clustermanager_u.SetHitManager(&hitmanager_u);
00484 
00485 //      Managed::HitManager hitmanager_v;
00486 //      Managed::ClusterManager clustermanager_v;
00487 //      clustermanager_v.SetHitManager(&hitmanager_v);
00488 
00489 
00490         MSG("ParticleFinder",Msg::kDebug) << "-event " << ievt << "  with " << evts[ievt]->nstrip << " strips" <<endl; 
00491 
00492         double strip_e=0;
00493 
00494 
00495         //Setup and give context
00496         DbiResultPtr<CalMIPCalibration> fResPtr;
00497         fResPtr.NewQuery(str->GetHeader().GetVldContext(),10);
00498 
00499         //Get Scale
00500         const CalMIPCalibration* mipcal = fResPtr.GetRow(0);
00501         double sigcorpermip= mipcal->GetScale();
00502 
00503         if(recoCnt<1)printf("Using SigCorPerMip %f\n",sigcorpermip);
00504 
00505 /*
00506         double sigcorpermip=0;
00507         if(n->GetHeader().GetVldContext().GetDetector()==Detector::kFar)sigcorpermip=kb->Get("FarSigCor/MIP");
00508         if(n->GetHeader().GetVldContext().GetDetector()==Detector::kNear)sigcorpermip=kb->Get("NearSigCor/MIP");
00509 */      
00510         //double pecutlevel=kb->Get("PECutLevel");
00511         //double MIPperPE=kb->Get("MIPperPE");
00512 
00513 
00514         //store the frontmost n hits in each view for containment
00515 
00516         double front_t[2][nfront];
00517         double front_z[2][nfront];
00518 
00519         for(int v=0;v<2;v++)
00520         for(int i=0;i<nfront;i++)
00521         {
00522                 front_t[v][i]=0;
00523                 front_z[v][i]=10000;
00524         }
00525 
00526         //Load processor with event strip data
00527         for(int i=0;i<evts[ievt]->nstrip;i++)
00528         {
00529         
00530                 double strip_z = stp[evts[ievt]->stp[i]]->z;
00531                 double strip_t = stp[evts[ievt]->stp[i]]->tpos;
00532                 int view = stp[evts[ievt]->stp[i]]->planeview;
00533         
00534                 if(! (view==2 || view ==3))continue;//we don't know what to do with any other view
00535                 
00536                 int forwardpos=nfront;
00537                 for(;forwardpos>0;forwardpos--)
00538                         if(front_z[view-2][forwardpos-1]<strip_z)break;
00539                         
00540                 //shift positions
00541                 for(int k=nfront-1;k>=forwardpos;k--)
00542                 {
00543                         front_z[view-2][k]=front_z[view-2][k-1];
00544                         front_t[view-2][k]=front_t[view-2][k-1];
00545                 }
00546                 
00547                 if(forwardpos<nfront)
00548                 {
00549                         front_z[view-2][forwardpos]=strip_z;
00550                         front_t[view-2][forwardpos]=strip_t;
00551                 }
00552                 
00553                         
00554         
00555                 //pe cut????
00556                 if((stp[evts[ievt]->stp[i]]->ph0.sigcor+stp[evts[ievt]->stp[i]]->ph1.sigcor)/sigcorpermip<pecutlevel*MIPperPE)continue;
00557         
00558 
00559                 int plane = stp[evts[ievt]->stp[i]]->plane;
00560                 int strip = stp[evts[ievt]->stp[i]]->strip;
00561                 double energy = (stp[evts[ievt]->stp[i]]->ph0.sigcor+stp[evts[ievt]->stp[i]]->ph1.sigcor)/sigcorpermip;
00562                 double tpos = stp[evts[ievt]->stp[i]]->tpos;
00563                 double z = stp[evts[ievt]->stp[i]]->z;
00564                 int planeview = stp[evts[ievt]->stp[i]]->planeview;
00565         
00566                 finder.AddStrip(plane, strip, energy, tpos, z, planeview);
00567                 strip_e+=energy;
00568 
00569 
00570                 hitmanager_u.InsertHit(planeview, plane, strip, z, tpos, energy);
00571 //              printf("adding strip %d %d %f\n",plane,strip,energy);
00572         }
00573 
00574 //      printf("stripenergy %f\n",strip_e);
00575         /*
00576         printf("front strips\n");
00577         for(int v=0;v<2;v++)
00578         {
00579                 printf("view %d\n",v+2);
00580         for(int i=0;i<nfront;i++)
00581         {
00582                 if(front_z[v][i]<10000)printf("(z %f t %f)",front_z[v][i],front_t[v][i]);
00583         }
00584                 printf("\n");
00585         }
00586         */
00587         
00588         CheckContainment(n, front_z, front_t);  
00589         //printf("forward contained: %d\n",n->event.contained);
00590 
00591         
00592 
00593 
00594                 
00595                 
00596 
00597 
00598 
00599         //clustermanager_u.DumpClusters();
00600 //      clustermanager_v.MakeClusters(0.2,0.05,0.01);
00601         //clustermanager_v.DumpClusters();
00602 
00603 
00604         
00605         Managed::ClusterSaver * cs=&n->cluster_saver;
00606 
00607         
00608         clustermanager_u.SetClusterSaver(cs);
00609 //      clustermanager_v.SetClusterSaver(cs);
00610 
00611         finder.SetClusterManagerU(clustermanager_u);
00612         finder.SetClusterManagerV(clustermanager_u);
00613         
00614 n->event.visenergy=strip_e;
00615 
00616         MSG("ParticleFinder",Msg::kDebug) << finder.GetStrips() <<" strips added"<<endl; 
00617 
00618         //process event, returns particleobject 
00619 
00620         
00621                 
00622 
00623         
00624         
00625         finder.Process(*n);
00626 
00627 
00628         
00629         //printf("p3d size %d\n",n->particles3d1->GetEntries());
00630 
00631         n->event.vtx_u=finder.vtx_u;
00632         n->event.vtx_v=finder.vtx_v;
00633         n->event.vtx_z=finder.vtx_z;    
00634         
00635 
00636         double x=(n->event.vtx_u-n->event.vtx_v)/sqrt(2.0);
00637         double y=(n->event.vtx_u+n->event.vtx_v)/sqrt(2.0);
00638         double z=n->event.vtx_z;
00639         
00640         //check for fiducial
00641 
00642         if(n->GetHeader().GetVldContext().GetDetector() ==Detector::kNear)
00643                 if(IsInsideNearFiducial_Nue_Standard(x,y,z,ismc)==1) 
00644                         n->event.inFiducial=1;  
00645  
00646         if(n->GetHeader().GetVldContext().GetDetector()==Detector::kFar)        
00647                 if(IsInsideFarFiducial_Nue_Standard(x,y,z,ismc)==1)
00648                         n->event.inFiducial=1;
00649 
00650 
00651 
00652 
00653 
00654                 n->SetName(outObjectName.c_str());
00655                 mom->AdoptFragment(n);
00656 
00657         recoCnt++;
00658 
00659 }
00660 
00661 
00662 
00663 
00664 
00666 //a check to see whats there...
00667 /*
00668  hft =( mom->GetFragmentList("ParticleObjectHolder"));
00669 n=dynamic_cast<ParticleObjectHolder *>(hft[0]);
00670 
00671 printf("#poh is %d   , #p3d %d\n",hft.size(),n->particles3d1->GetEntries());
00672 */
00673 
00674 
00675 
00677 
00678 
00679 
00680   return JobCResult::kPassed; // kNoDecision, kFailed, etc.
00681 }

void ParticleFinder::Reset (  )  [virtual]

Implement to reset oneself

(Document me!)

Reimplemented from JobCModule.

Definition at line 846 of file ParticleFinder.cxx.

00847 {
00851 }


Member Data Documentation

Definition at line 80 of file ParticleFinder.h.

Referenced by Reco().

int ParticleFinder::DoMRCC [private]

Definition at line 84 of file ParticleFinder.h.

Referenced by Config(), ParticleFinder(), and Reco().

Definition at line 59 of file ParticleFinder.h.

Referenced by Reco().

std::string ParticleFinder::kPOTTreeName [private]

Definition at line 69 of file ParticleFinder.h.

Referenced by Reco().

int ParticleFinder::lastrun [private]

Definition at line 75 of file ParticleFinder.h.

Referenced by FillPot(), and ParticleFinder().

double ParticleFinder::MEUperGeV [private]

Definition at line 88 of file ParticleFinder.h.

Referenced by Config(), and Reco().

double ParticleFinder::MIPperPE [private]

Definition at line 87 of file ParticleFinder.h.

Referenced by Config(), ParticleFinder(), and Reco().

TF1 * ParticleFinder::osceq = 0 [static, private]

Definition at line 64 of file ParticleFinder.h.

Referenced by ParticleFinder(), Reco(), and ~ParticleFinder().

std::string ParticleFinder::outObjectName [private]

Definition at line 83 of file ParticleFinder.h.

Referenced by Config(), and Reco().

double ParticleFinder::pecutlevel [private]

Definition at line 86 of file ParticleFinder.h.

Referenced by Config(), ParticleFinder(), and Reco().

int ParticleFinder::recoCnt [private]

Definition at line 81 of file ParticleFinder.h.

Referenced by ParticleFinder(), and Reco().

Definition at line 71 of file ParticleFinder.h.

Referenced by ParticleFinder(), Reco(), and ~ParticleFinder().


The documentation for this class was generated from the following files:

Generated on 22 Nov 2017 for loon by  doxygen 1.6.1