CondensedNtpModuleAtm Class Reference

#include <CondensedNtpModuleAtm.h>

Inheritance diagram for CondensedNtpModuleAtm:

List of all members.

Public Member Functions
	CondensedNtpModuleAtm ()
virtual	~CondensedNtpModuleAtm ()
void	BeginJob ()
JobCResult	Ana (const MomNavigator *mom)
const Registry &	DefaultConfig () const
void	Config (const Registry &r)
void	Help ()
void	EndJob ()
Private Member Functions
void	FillEventInformation (ANtpRecoNtpManipulator *ntpManipulator, NtpSREvent *ntpEvent)
void	FillMCInformation (NtpMCTruth *ntpMCTruth, NtpMCStdHep *ntpMCStdHep)
void	FillTrackInformation (NtpSRTrack *ntpTrack, TClonesArray *stripArray, double azimuth, double ra, double dec)
void	FillTrackTimingVariables (NtpSRTrack *ntpTrack, TClonesArray *stripArray)
void	FillTrackMagneticBendingVariables (NtpSRTrack *ntpTrack, TClonesArray *stripArray)
void	ResetTreeVariables ()
Double_t	GetTimeWeight (Float_t ph)
Double_t	ArrivalTime_Weight (Double_t npe) const
Double_t	ArrivalTime_Uncertainty (Double_t npe) const
void	WeightedAverage (const Int_t n, const Double_t *x, const Double_t *w, Double_t *parm, Double_t *eparm)
void	LinearFit_Weighted (const Int_t n, const Double_t *x, const Double_t *y, const Double_t *w, Double_t *parm, Double_t *eparm)
void	FitMinimizingResidual (Int_t &n, Double_t *x, Double_t *y, Double_t *w, Double_t *param, Double_t *eparam, bool findSlope=true)
void	FindLinearFit (Double_t *x, Double_t *y, Double_t *weight, Int_t nPoints, Double_t &a1, Double_t &a2, Double_t &chiSq)
void	FindChi2 (Double_t *x, Double_t *y, Double_t *weight, Int_t nPoints, Double_t &a1, Double_t &a2, Double_t &chiSq)
void	FindStraightLineDeviation (Double_t *x, Double_t *y, Int_t nPoints, Double_t &a1, Double_t &a2, Double_t &deviation)
void	FindMeanAndRMS (double *x, double *weight, int nPoints, double &mean, double &rms)
Float_t	CheapBaseline (float cosz)
Private Attributes
std::string	fFileName
std::string	fRockMapFileName
std::string	fTreeName
std::string	fBeamTreeName
std::string	fBeamPath
TFile *	fNtpFile
TTree *	fNtuple
TTree *	fFailTree
TTree *	fBeamTree
Int_t	fDetector
Int_t	fGoodTrack
Int_t	fDataType
Int_t	fEndPlane
Double_t	fHornCurrent
Double_t	fCurrentPOT
Double_t	fBeamHPos
Double_t	fBeamVPos
Double_t	fTargetPos
Double_t	fBeamTimeSec
Int_t	fIsMultiple
Int_t	fValidPlaneFail
Int_t	fFailDeMux
Int_t	fMajorRelease
Double_t	fPathCosZenDeg [1500]
Double_t	fPathAzimuthDeg [1500]
Double_t	fColDenCosZen [1500]
Double_t	fDensity [1500]
bool	fNewMC
ANtpEventInfo *	fEventInfo
ANtpHeaderInfo *	fHeaderInfo
ANtpBeamInfo *	fBeamInfo
ANtpTrackInfoAtm *	fTrackInfo
ANtpTruthInfoAtm *	fTruthInfo
ANtpInfoObjectFillerBeam *	fInfoFiller
Int_t	fCtr

---

Constructor & Destructor Documentation

CondensedNtpModuleAtm::CondensedNtpModuleAtm (   )

Definition at line 70 of file CondensedNtpModuleAtm.cxx. References MSG. : fFileName("analysisNtuple.root"), fTreeName("analysisNtuple"), fNtpFile(0), fNtuple(0), fBeamTree(0), fDetector(1), fEndPlane(485), fEventInfo(0), fHeaderInfo(0), fTrackInfo(0), fTruthInfo(0), fCtr(0) { MSG("JobC", Msg::kDebug) << "CondensedNtpModuleAtm::Constructor" << endl; fHeaderInfo = new ANtpHeaderInfo(); fBeamInfo = new ANtpBeamInfo(); fEventInfo = new ANtpEventInfo(); fTrackInfo = new ANtpTrackInfoAtm(); fTruthInfo = new ANtpTruthInfoAtm(); fInfoFiller = new ANtpInfoObjectFillerBeam(); }

CondensedNtpModuleAtm::~CondensedNtpModuleAtm (   )  [virtual]

Definition at line 97 of file CondensedNtpModuleAtm.cxx. References MSG. { MSG("JobC", Msg::kDebug) << "CondensedNtpModuleAtm::Destructor" << endl; if(fHeaderInfo) delete fHeaderInfo; if(fBeamInfo) delete fBeamInfo; if(fEventInfo) delete fEventInfo; if(fTrackInfo) delete fTrackInfo; if(fTruthInfo) delete fTruthInfo; }

---

Member Function Documentation

JobCResult CondensedNtpModuleAtm::Ana (  const MomNavigator *  mom  )  [virtual]

Implement this for read only access to the MomNavigator Reimplemented from JobCModule. Definition at line 190 of file CondensedNtpModuleAtm.cxx. References NtpSRCosmicRay::azimuth, NtpSRCosmicRay::dec, ANtpHeaderInfo::detector, fBeamInfo, fDetector, fEndPlane, fEventInfo, fFailDeMux, fFailTree, fGoodTrack, fHeaderInfo, ANtpInfoObjectFillerBeam::FillBeamInformation(), FillEventInformation(), ANtpInfoObjectFiller::FillHeaderInformation(), FillMCInformation(), FillTrackInformation(), FillTrackMagneticBendingVariables(), FillTrackTimingVariables(), fInfoFiller, fIsMultiple, fMajorRelease, fNtuple, fTrackInfo, fValidPlaneFail, ANtpRecoNtpManipulator::GetCosmicRayInfo(), ANtpRecoNtpManipulator::GetDmxStatusInfo(), ANtpEventManipulator::GetEvent(), ANtpRecoNtpManipulator::GetEventManipulator(), MomNavigator::GetFragment(), ANtpRecoNtpManipulator::GetMCManipulator(), ANtpMCManipulator::GetNtpMCStdHep(), ANtpMCManipulator::GetNtpMCTruth(), ANtpMCManipulator::GetNtpTHEvent(), ANtpEventManipulator::GetPrimaryTrack(), ANtpRecoNtpManipulator::GetSnarlEventSummary(), ANtpRecoNtpManipulator::GetStripArray(), VldContext::GetTimeStamp(), RecRecordImp< T >::GetVldContext(), ANtpEventInfo::index, NtpSRDmxStatus::ismultimuon, MSG, NtpTHEvent::neumc, NtpSREventSummary::nevent, NtpSRDmxStatus::nonphysicalfail, NtpSREventSummary::ntrack, ANtpTrackInfo::planes, NtpSRCosmicRay::ra, ResetTreeVariables(), ANtpHeaderInfo::run, ANtpEventManipulator::SetEventInSnarl(), JobCResult::SetFailed(), ANtpRecoNtpManipulator::SetPrimaryShowerCriteria(), ANtpRecoNtpManipulator::SetPrimaryTrackCriteria(), ANtpInfoObjectFiller::SetStripArray(), JobCResult::SetWarning(), ANtpHeaderInfo::snarl, ANtpTrackInfo::totalStrips, ANtpEventInfo::totalStrips, NtpSRDmxStatus::ustrayplanes, NtpSRDmxStatus::uvalidplanes, NtpSRDmxStatus::validplanesfail, NtpSRDmxStatus::vertexplanefail, NtpSRDmxStatus::vstrayplanes, NtpSRDmxStatus::vvalidplanes, and ANtpHeaderInfo::year. { MSG("CondensedNtpModuleAtm", Msg::kDebug) << "start ana" << endl; JobCResult result(JobCResult::kPassed); //get the records from MOM assert(mom); NtpSRRecord *record = dynamic_cast<NtpSRRecord *>(mom->GetFragment("NtpSRRecord")); NtpStRecord *stRecord = dynamic_cast<NtpStRecord *>(mom->GetFragment("NtpStRecord")); if(!record && !stRecord){ MSG("CondensedNtpModuleAtm", Msg::kWarning) << "Could not get NtpSR or NtpSt Record from MOM" << endl; result.SetWarning().SetFailed(); return result; }//end if no NtpSRRecord NtpMCRecord *mcRecord = 0; NtpTHRecord *thRecord = 0; mcRecord = dynamic_cast<NtpMCRecord *>(mom->GetFragment("NtpMCRecord")); thRecord = dynamic_cast<NtpTHRecord *>(mom->GetFragment("NtpTHRecord")); //get the record for the beam tree NtpBDLiteRecord *bdRecord = dynamic_cast<NtpBDLiteRecord *>(mom->GetFragment("NtpBDLiteRecord")); if(!bdRecord) MSG("CondensedNtpModuleAtm", Msg::kDebug) << "NO BEAM RECORD" << endl; //fill the header information DetectorType::Detector_t det = DetectorType::kFar; if(fDetector<1) det = DetectorType::kNear; //get a sim flag SimFlag::SimFlag_t dataType = SimFlag::kData; if(fDataType) dataType = SimFlag::kMC; //instantiate a NtpHelper object to help you get the info you want ANtpRecoNtpManipulator *ntpManipulator = 0; if(record) ntpManipulator = new ANtpRecoNtpManipulator(record, mcRecord, thRecord); else if(stRecord) ntpManipulator = new ANtpRecoNtpManipulator(stRecord); VldTimeStamp timeStamp = stRecord ? stRecord->GetVldContext()->GetTimeStamp() : record->GetVldContext()->GetTimeStamp(); MSG("CondensedNtpModuleAtm", Msg::kDebug) << timeStamp << endl; // TStopwatch stopwatch; // stopwatch.Start(); fInfoFiller->SetStripArray(ntpManipulator->GetStripArray()); fInfoFiller->FillHeaderInformation(ntpManipulator, det, dataType, fHeaderInfo); // stopwatch.Stop(); // MSG("CondensedNtpModuleAtm", Msg::kDebug) << "took " << stopwatch.RealTime() // << " to fill header" << endl; // stopwatch.Reset(); // stopwatch.Start(); if(fHeaderInfo->year > 2004){ if(bdRecord) fInfoFiller->FillBeamInformation(bdRecord, fBeamInfo); else fInfoFiller->FillBeamInformation(timeStamp, det, dataType, fBeamInfo); } // stopwatch.Stop(); // MSG("CondensedNtpModuleAtm", Msg::kDebug) << "took " << stopwatch.RealTime() // << " to fill beam" << endl; // stopwatch.Reset(); //set up which flags you want to use to determine the primary shower or track //a value of 0 for a flag means it will not be used ntpManipulator->SetPrimaryTrackCriteria(0,1,0); // nplanes, length, total pulse height ntpManipulator->SetPrimaryShowerCriteria(0,1); // nplanes, total pulse height //check that the current event passed demuxing bool passedDeMux = true; fGoodTrack = true; fIsMultiple = 0; fValidPlaneFail = 0; fFailDeMux = 0; MSG("CondensedNtpModuleAtm", Msg::kDebug) << "run,snarl = " << fHeaderInfo->run << "," << fHeaderInfo->snarl << " " << ntpManipulator->GetSnarlEventSummary().ntrack << " " << passedDeMux << endl; if(fHeaderInfo->detector == DetectorType::kFar){ if(ntpManipulator->GetDmxStatusInfo().ismultimuon){ passedDeMux = false; fIsMultiple = 1; } if(ntpManipulator->GetDmxStatusInfo().validplanesfail){ passedDeMux = false; fValidPlaneFail = 1; } if(fMajorRelease < 18){ if(ntpManipulator->GetDmxStatusInfo().vertexplanefail || ntpManipulator->GetDmxStatusInfo().nonphysicalfail) passedDeMux = false; UShort_t uValidPlanes = ntpManipulator->GetDmxStatusInfo().uvalidplanes; UShort_t uStrayPlanes = ntpManipulator->GetDmxStatusInfo().ustrayplanes; UShort_t vValidPlanes = ntpManipulator->GetDmxStatusInfo().vvalidplanes; UShort_t vStrayPlanes = ntpManipulator->GetDmxStatusInfo().vstrayplanes; Float_t ufrac = 0., vfrac = 0.; if(uValidPlanes>0) ufrac = (1.*uStrayPlanes)/(1.*uValidPlanes); else ufrac = 1.; if(vValidPlanes>0) vfrac = (1.*vStrayPlanes)/(1.*vValidPlanes); else vfrac = 1.; if(ufrac>=0.45) passedDeMux = false; else if( vfrac>=0.45 ) passedDeMux = false; else if(uStrayPlanes>=3&&uValidPlanes<20) passedDeMux = false; else if(vStrayPlanes>=3&&vValidPlanes<20) passedDeMux = false; else if(uStrayPlanes>=4&&uValidPlanes>=20 && uValidPlanes<=40) passedDeMux = false; else if(vStrayPlanes>=4&&vValidPlanes>=20 && vValidPlanes<=40) passedDeMux = false; else if(ufrac>=0.1 && uValidPlanes>=40) passedDeMux = false; else if(vfrac>=0.1 && vValidPlanes>=40) passedDeMux = false; }//end if before switch over to cambridge demuxing }//end if far detector //if the event isnt demuxed well, then forget it also forget it //if there arent any tracks in the event if(!passedDeMux || ntpManipulator->GetSnarlEventSummary().ntrack==0 || fEndPlane<248){ fFailTree->Fill(); return result; } //set up some pointers to typical NtpSR objects NtpSRTrack *ntpTrack = 0; NtpSREvent *ntpEvent = 0; NtpMCTruth *ntpMCTruth = 0; NtpMCStdHep *ntpMCStdHep = 0; NtpTHEvent *ntpTHEvent = 0; //loop over the events in the snarl for(Int_t i = 0; i < ntpManipulator->GetSnarlEventSummary().nevent; ++i){ MSG("CondensedNtpModuleAtm", Msg::kDebug) << "on event " << i << endl; ResetTreeVariables(); fEventInfo->index = i; // stopwatch.Start(); //get event i. the call to ANtpRecoNtpManipulator::GetEventInSnarl sets the //NtpSREvent data member in that object to the current event so that //you can later call for the primary track, shower, etc. ntpManipulator->GetEventManipulator()->SetEventInSnarl(i); ntpEvent = ntpManipulator->GetEventManipulator()->GetEvent(); if(ntpEvent) FillEventInformation(ntpManipulator,ntpEvent); // stopwatch.Stop(); // MSG("CondensedNtpModuleAtm", Msg::kInfo) << "took " << stopwatch.RealTime() // << " to fill event" << endl; // stopwatch.Reset(); // stopwatch.Start(); //get the primary track for the event - if no track is present it //returns 0 fGoodTrack = false; ntpTrack = ntpManipulator->GetEventManipulator()->GetPrimaryTrack(); if(ntpTrack){ fGoodTrack = true; FillTrackInformation(ntpTrack, ntpManipulator->GetStripArray(), ntpManipulator->GetCosmicRayInfo().azimuth, ntpManipulator->GetCosmicRayInfo().ra, ntpManipulator->GetCosmicRayInfo().dec); FillTrackMagneticBendingVariables(ntpTrack, ntpManipulator->GetStripArray()); if(fDetector>0&&fGoodTrack) FillTrackTimingVariables(ntpTrack, ntpManipulator->GetStripArray()); } // stopwatch.Stop(); // MSG("CondensedNtpModuleAtm", Msg::kInfo) << "took " << stopwatch.RealTime() // << " to fill track" << endl; // stopwatch.Reset(); // stopwatch.Start(); // Get best neu match from Truthhelper ntpTHEvent = ntpManipulator->GetMCManipulator()->GetNtpTHEvent(i); MSG("CondensedNtpModuleAtm", Msg::kDebug) << "start filling mc info " << ntpTHEvent->neumc << endl; if(ntpTHEvent) ntpMCTruth = ntpManipulator->GetMCManipulator()->GetNtpMCTruth(ntpTHEvent->neumc); if(ntpMCTruth){ //get the NtpMCStdHep object - use the fact that i put the values i want in the parent //place, which means i want the 0th entry if(fNewMC) ntpMCStdHep = ntpManipulator->GetMCManipulator()->GetNtpMCStdHep(0); FillMCInformation(ntpMCTruth, ntpMCStdHep); } // stopwatch.Stop(); // MSG("CondensedNtpModuleAtm", Msg::kInfo) << "took " << stopwatch.RealTime() // << " to fill truth" << endl; // stopwatch.Reset(); if(fEventInfo->totalStrips<fTrackInfo->totalStrips) fGoodTrack = false; if(fTrackInfo->planes<1) fGoodTrack = false; if(fTrackInfo->totalStrips<1) fGoodTrack = false; MSG("CondensedNtpModuleAtm", Msg::kDebug) << "strips = " << fEventInfo->totalStrips << "/" << fTrackInfo->totalStrips << " planes = " << fTrackInfo->planes << " good track = " << (int)fGoodTrack << endl; //got all the necessary sr info, now fill the ntuple if(!fGoodTrack) continue; fNtuple->Fill(); }//end loop over events for this snarl ++fCtr; delete ntpManipulator; return result; }

Double_t CondensedNtpModuleAtm::ArrivalTime_Uncertainty (  Double_t  npe  )  const [private]

Definition at line 1025 of file CondensedNtpModuleAtm.cxx. Referenced by ArrivalTime_Weight(). { // cout << "in ArrivalTime_Uncertainty" << endl; Double_t sigma = 0.; // sigma is uncertainty in arrival of first p.e. // a simple monte carlo was done to calculate the uncertainty in the // arrival time of the first p.e. as a function of # of p.e. // // the only assumption made was about the shape of the pulse, a 2 ns risetime // and an 8 ns decay time of the fluor // // rlee feb 2002 // if(npe<10.) sigma = 0.90+8.0/npe; else sigma = 5.0*exp(-0.5*log(npe)); // cout << sigma << endl; return sigma; }

Double_t CondensedNtpModuleAtm::ArrivalTime_Weight (  Double_t  npe  )  const [private]

Definition at line 1013 of file CondensedNtpModuleAtm.cxx. References ArrivalTime_Uncertainty(). Referenced by GetTimeWeight(). { // cout << "in ArrivalTime_Weight" << endl; Double_t weight = 1./ArrivalTime_Uncertainty(npe); // //cout << weight << endl; return weight*weight; }

void CondensedNtpModuleAtm::BeginJob (   )  [virtual]

Implement for notification of begin of job Reimplemented from JobCModule. Definition at line 111 of file CondensedNtpModuleAtm.cxx. References fBeamInfo, fBeamPath, fBeamTreeName, fColDenCosZen, fDataType, fDensity, fEventInfo, fFailDeMux, fFailTree, fFileName, fHeaderInfo, fIsMultiple, fNtpFile, fNtuple, fPathAzimuthDeg, fPathCosZenDeg, fRockMapFileName, fTrackInfo, fTreeName, fTruthInfo, fValidPlaneFail, and MSG. { // save the current working directory TDirectory* savedir = gDirectory; //create the new TFile for holding the electronics tree // this changes the value of gDirectory fNtpFile = new TFile(fFileName.c_str(),"RECREATE"); MSG("CondensedNtpModuleAtm", Msg::kDebug) << fFileName << " " << fTreeName << " " << fNtpFile->GetName() << endl; // create TTree, these will attach themselves to the current //working directory fNtuple = new TTree(fTreeName.c_str(), "Analysis Tree"); fNtuple->Branch("header.", "ANtpHeaderInfo", &fHeaderInfo, 64000, 2); fNtuple->Branch("beam.", "ANtpBeamInfo", &fBeamInfo, 64000, 2); fNtuple->Branch("event.", "ANtpEventInfo", &fEventInfo, 64000, 2); fNtuple->Branch("track.", "ANtpTrackInfoAtm", &fTrackInfo, 64000, 2); if(fDataType==1) fNtuple->Branch("truth.", "ANtpTruthInfoAtm", &fTruthInfo, 64000, 2); MSG("CondensedNtpModuleAtm", Msg::kDebug) << "got branches" << endl; fFailTree = new TTree("failures", "Failed Events"); fFailTree->Branch("header.", "ANtpHeaderInfo", &fHeaderInfo, 64000, 2); fFailTree->Branch("multiple", &fIsMultiple, "multiple/I"); fFailTree->Branch("validplane", &fValidPlaneFail, "validplane/I"); fFailTree->Branch("failDeMux", &fFailDeMux, "failDeMux/I"); // change back to current working directory before leaving constructor savedir->cd(); //load the vertical intensity information into some arrays //get the slant depths and the number of events for each slant depth //read in rock data // http://www.hep.umn.ed/~schubert/far/s2rock/slantdepth.data // 1-cos(zen); azi(degrees); slantdepth (gmcm**2) * cos(zen); // density (gm/cm**2) // N.B., units for density are incorrectly labelled in the file ifstream densityFile(fRockMapFileName.c_str()); if( !densityFile.is_open() ){ cout << "could not open rock map - exit" << endl; exit(0); } int ir = 0; while( ir<1404 ){ densityFile >> fPathCosZenDeg[ir] >> fPathAzimuthDeg[ir] >> fColDenCosZen[ir] >> fDensity[ir]; ++ir; } //get the beam info tree //check to see if you need to calculate the total Protons on Target MSG("BeamTestModule", Msg::kInfo) << "beam info path = " << fBeamPath.c_str() << " tree name = " << fBeamTreeName.c_str() << endl; // if(fDataType < 1){ // TFile *beamFile = new TFile(fBeamPath.c_str()); // fBeamTree = dynamic_cast<TTree *>(beamFile->Get(fBeamTreeName.c_str())); // //build the index for the tree to select the utc timestamps // fBeamTree->BuildIndex("utc"); // fBeamTree->SetBranchAddress("intensity", &fCurrentPOT); // fBeamTree->SetBranchAddress("hornCurrent", &fHornCurrent); // fBeamTree->SetBranchAddress("hpos", &fBeamHPos); // fBeamTree->SetBranchAddress("vpos", &fBeamVPos); // fBeamTree->SetBranchAddress("tgtpos", &fTargetPos); // } return; }

Float_t CondensedNtpModuleAtm::CheapBaseline (  float  cosz  )  [private]

Definition at line 1399 of file CondensedNtpModuleAtm.cxx. Referenced by FillMCInformation(). { static float Rearth = 6378140.0; /* Equatorial radius of earth */ static float nuHeight = 1500.0; /* Average Neutrino production height */ static float altitudeSK = -210.0; /* Altitude of the detector - MINOS is 689ft below sea level*/ static float x1s = (Rearth+nuHeight)*(Rearth+nuHeight); static float x2 = Rearth+altitudeSK; static float x2s = (Rearth+altitudeSK)*(Rearth+altitudeSK); return sqrt(x1s+x2s*(cosz*cosz-1.0))-x2*cosz; }

void CondensedNtpModuleAtm::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. Reimplemented from JobCModule. Definition at line 474 of file CondensedNtpModuleAtm.cxx. References fBeamPath, fBeamTreeName, fDataType, fDetector, fEndPlane, fFileName, fMajorRelease, fNewMC, fRockMapFileName, fTreeName, and Registry::Get(). { int tmpb; // a temp bool. See comment under DefaultConfig... int tmpi; // a temp int. double tmpd; // a temp double. const char* tmps; // a temp string if (r.Get("FileName", tmps)) fFileName = tmps; if (r.Get("RockMapFileName", tmps)) fRockMapFileName = tmps; if (r.Get("TreeName", tmps)) fTreeName = tmps; if (r.Get("BeamPath", tmps)) fBeamPath = tmps; if (r.Get("BeamTreeName", tmps)) fBeamTreeName = tmps; if (r.Get("EndPlane", tmpi)) fEndPlane = tmpi; if (r.Get("Detector", tmpi)) fDetector = tmpi; if (r.Get("DataType", tmpi)) fDataType = tmpi; if (r.Get("NewMC", tmpb)) fNewMC = tmpb; if (r.Get("MajorRelease", tmpi)) fMajorRelease = tmpi; //quiet the compiler warnings tmpb = 0; tmpd = 1.; return; }

const Registry & CondensedNtpModuleAtm::DefaultConfig (   )  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; } Reimplemented from JobCModule. Definition at line 446 of file CondensedNtpModuleAtm.cxx. References Registry::LockValues(), Registry::Set(), and Registry::UnLockValues(). { int itrue = 1; // work around for lack of bool in registry int ifalse = 0; // work around for lack of bool in registry static Registry r; r.UnLockValues(); r.Set("FileName", "analysisTree.root"); r.Set("RockMapFileName", "density.txt"); r.Set("TreeName", "analysisTree"); r.Set("BeamPath", "/afs/fnal.gov/files/data/minos/d04/brebel/beamInfo.root"); r.Set("BeamTreeName", "beam"); r.Set("EndPlane", 485); r.Set("Detector", 1); r.Set("DataType", 0); r.Set("NewMC", itrue); r.Set("MajorRelease", 14); r.LockValues(); itrue = ifalse; //quiet the compiler warnings return r; }

void CondensedNtpModuleAtm::EndJob (   )  [virtual]

Implement for notification of end of job Reimplemented from JobCModule. Definition at line 420 of file CondensedNtpModuleAtm.cxx. References fFailTree, fNtpFile, fNtuple, and MSG. { MSG("CondensedNtpModuleAtm", Msg::kDebug) << "start end job method" << endl; //save current working directory TDirectory *savedir = gDirectory; //cd to the "directory" of the file fNtpFile->cd(); MSG("CondensedNtpModuleAtm", Msg::kInfo) << fNtuple->GetEntries() << endl; fNtuple->Write(); fFailTree->Write(); //go back to original working directory savedir->cd(); fNtpFile->Write(); fNtpFile->Close(); return; }

void CondensedNtpModuleAtm::FillEventInformation (  ANtpRecoNtpManipulator *  ntpManipulator, NtpSREvent *  ntpEvent )  [private]

Definition at line 935 of file CondensedNtpModuleAtm.cxx. References fEventInfo, ANtpInfoObjectFiller::FillEventInformation(), fInfoFiller, and MSG. Referenced by Ana(). { MSG("CondensedNtpModuleAtm", Msg::kDebug) << "in FillEventInformation" << endl; fInfoFiller->FillEventInformation(ntpManipulator, ntpEvent, fEventInfo); return; }

void CondensedNtpModuleAtm::FillMCInformation (  NtpMCTruth *  ntpMCTruth, NtpMCStdHep *  ntpMCStdHep )  [private]

Definition at line 949 of file CondensedNtpModuleAtm.cxx. References ANtpTruthInfoAtm::azimuth, ANtpTruthInfoAtm::baseLine, CheapBaseline(), ANtpInfoObjectFiller::FillMCTruthInformation(), fInfoFiller, fTruthInfo, ANtpTruthInfo::leptonDCosX, ANtpTruthInfo::leptonDCosY, ANtpTruthInfo::leptonDCosZ, MSG, ANtpTruthInfo::nuEnergy, NtpMCTruth::p4neu, NtpMCStdHep::vtx, and ANtpTruthInfoAtm::weight. Referenced by Ana(). { MSG("CondensedNtpModuleAtm", Msg::kDebug) << "in FillMCInformation" << endl; const double rotMatrix[] = { 0.894507011,0,0.447053919, 0, 1, 0, -0.447053919,0,0.894507011}; double dcosx_ideal, dcosy_ideal, dcosz_ideal; fInfoFiller->FillMCTruthInformation(ntpMCTruth, fTruthInfo); if(fTruthInfo->nuEnergy>0.) fTruthInfo->baseLine = 0.001*CheapBaseline(-ntpMCTruth->p4neu[1] /TMath::Abs(fTruthInfo->nuEnergy)); //get the azimuth from AstroUtill dcosx_ideal = (rotMatrix[0]*fTruthInfo->leptonDCosX + rotMatrix[1]*fTruthInfo->leptonDCosY + rotMatrix[2]*fTruthInfo->leptonDCosZ); dcosy_ideal = (rotMatrix[3]*fTruthInfo->leptonDCosX + rotMatrix[4]*fTruthInfo->leptonDCosY + rotMatrix[5]*fTruthInfo->leptonDCosZ); dcosz_ideal = (rotMatrix[6]*fTruthInfo->leptonDCosX + rotMatrix[7]*fTruthInfo->leptonDCosY + rotMatrix[8]*fTruthInfo->leptonDCosZ); fTruthInfo->azimuth = 180.*TMath::ATan2(-dcosx_ideal,dcosz_ideal)/TMath::Pi(); if(fTruthInfo->azimuth < 0.) fTruthInfo->azimuth += 360.; //use the NtpMCStdHep object to figure out the weight for this event based on the //fact that the new MC picks events at the surface. we plugged the value into //the vtx Y location in the stdhep common bloc, gminos multiplies it by 10^-3 //because it thought it was converting mm to m if(ntpMCStdHep){ fTruthInfo->weight = 1.e3*ntpMCStdHep->vtx[1]; } return; }

void CondensedNtpModuleAtm::FillTrackInformation (  NtpSRTrack *  ntpTrack, TClonesArray *  stripArray, double  azimuth, double  ra, double  dec )  [private]

Definition at line 501 of file CondensedNtpModuleAtm.cxx. References ANtpTrackInfoAtm::alternateChi2, ANtpTrackInfoAtm::azimuth, NtpSRTrackTime::cdtds, NtpSRVertex::dcosx, ANtpTrackInfo::dcosXEnd, ANtpTrackInfo::dcosXVtx, NtpSRVertex::dcosy, ANtpTrackInfo::dcosYEnd, ANtpTrackInfo::dcosYVtx, NtpSRVertex::dcosz, ANtpTrackInfo::dcosZEnd, ANtpTrackInfo::dcosZVtx, ANtpTrackInfoAtm::dec, NtpSRFiducial::dr, NtpSRFiducial::dz, NtpSRTrack::end, ANtpTrackInfo::endX, ANtpTrackInfo::endY, ANtpTrackInfo::endZ, fDensity, fDetector, fEndPlane, fEventInfo, fGoodTrack, NtpSRTrack::fidend, ANtpTrackInfoAtm::fiducialEndDr, ANtpTrackInfoAtm::fiducialEndDz, ANtpTrackInfoAtm::fiducialVtxDr, ANtpTrackInfoAtm::fiducialVtxDz, NtpSRTrack::fidvtx, ANtpInfoObjectFiller::FillTrackInformation(), fInfoFiller, NtpSRTrack::fit, fTrackInfo, ANtpTrackInfoAtm::impactParameter, ANtpTrackInfoAtm::invBeta, ANtpTrackInfo::length, NtpSRTrack::momentum, MSG, NtpSRStrip::ndigit, NtpSRFitTrack::ndof, NtpSRTrack::nstrip, NtpSRTrackPlane::ntrklike, NtpSRPlane::nu, NtpSRPlane::nv, NtpSRStrip::plane, NtpSRTrack::plane, ANtpEventInfo::planes, ANtpTrackInfoAtm::planesIn10Meter, ANtpTrackInfoAtm::planesIn15Meter, ANtpTrackInfoAtm::planesIn20Meter, ANtpTrackInfoAtm::planesIn25Meter, ANtpTrackInfoAtm::planesIn30Meter, ANtpTrackInfoAtm::planesIn35Meter, ANtpTrackInfoAtm::planesIn40Meter, ANtpTrackInfoAtm::planeUseFraction, NtpSRStrip::planeview, ANtpTrackInfo::pulseHeight, ANtpEventInfo::pulseHeight, NtpSRMomentum::qp, ANtpTrackInfoAtm::ra, ANtpTrackInfoAtm::signalUseFraction, ANtpTrackInfoAtm::slantDepth, NtpSRTrack::stp, NtpSRTrack::stpu, NtpSRTrack::stpv, NtpSRTrack::stpx, NtpSRTrack::stpy, NtpSRStrip::strip, NtpSRTrack::time, ANtpTrackInfo::totalStrips, NtpSRStrip::tpos, ANtpTrackInfoAtm::trackLikePlanes, ANtpTrackInfoAtm::twoEndStripFraction, ANtpTrackInfoAtm::uvAsymmetry, NtpSRTrack::vtx, ANtpTrackInfo::vtxX, ANtpTrackInfo::vtxY, ANtpTrackInfo::vtxZ, NtpSRVertex::x, NtpSRVertex::y, and NtpSRVertex::z. Referenced by Ana(). { MSG("CondensedNtpModuleAtm", Msg::kDebug) << "filling track info" << endl; //fill the basics first fInfoFiller->FillTrackInformation(ntpTrack, fTrackInfo); NtpSRStrip *ntpStrip = 0; Int_t plane = 0, strip = 0, numDoubleEndedStrips = 0, prevPlane = -1; float xPos = 0., yPos = 0.; Float_t radius = 0.; Int_t index = 0; //stuff to find the slant depth double vec1[39],vec2[36]; int index1, index2, index3; for(int i=0;i<39;++i){ vec1[i] = i*0.02; // cout << vec1[i] << endl; if(i<36){ vec2[i] = i*10. + 5.; // cout << vec2[i] << endl; } } // Kashahara picks events according to where they come from // she uses altitude, not zenith angle index1 = TMath::BinarySearch(39,vec1,1.-(-fTrackInfo->dcosYVtx))+1; index2 = TMath::BinarySearch(36,vec2,azimuth)+1; index3 = (index1)*36 +(index2); fTrackInfo->slantDepth = fDensity[index3]*0.01; //density is in 10^5g/cm^2, ie km.w.e, //0.01 converts to m.w.e 10^3m/km/(10^5g/cm^2/km.w.e) fTrackInfo->azimuth = azimuth; fTrackInfo->ra = ra; fTrackInfo->dec = dec; fTrackInfo->trackLikePlanes = ntpTrack->plane.ntrklike; Int_t nVPlanes = ntpTrack->plane.nv; Int_t nUPlanes = ntpTrack->plane.nu; fTrackInfo->invBeta = ntpTrack->time.cdtds; fTrackInfo->planesIn10Meter = 0; fTrackInfo->planesIn15Meter = 0; fTrackInfo->planesIn20Meter = 0; fTrackInfo->planesIn25Meter = 0; fTrackInfo->planesIn30Meter = 0; fTrackInfo->planesIn35Meter = 0; fTrackInfo->planesIn40Meter = 0; //check if the zenith angle indicates the right direction for the near detector. //if not, swap the ends and adjust the azimuth if(fTrackInfo->dcosYVtx > 0. && fDetector<1){ fTrackInfo->vtxX = ntpTrack->end.x; fTrackInfo->vtxY = ntpTrack->end.y; fTrackInfo->vtxZ = ntpTrack->end.z; fTrackInfo->endX = ntpTrack->vtx.x; fTrackInfo->endY = ntpTrack->vtx.y; fTrackInfo->endZ = ntpTrack->vtx.z; fTrackInfo->dcosXVtx = -ntpTrack->end.dcosx; fTrackInfo->dcosYVtx = -ntpTrack->end.dcosy; fTrackInfo->dcosZVtx = -ntpTrack->end.dcosz; fTrackInfo->dcosXEnd = -ntpTrack->vtx.dcosx; fTrackInfo->dcosYEnd = -ntpTrack->vtx.dcosy; fTrackInfo->dcosZEnd = -ntpTrack->vtx.dcosz; const double rotMatrix[] = { 0.894507011,0,0.447053919, 0, 1, 0, -0.447053919,0,0.894507011}; double dcosx_ideal, dcosy_ideal, dcosz_ideal; //get the azimuth from AstroUtill dcosx_ideal = (rotMatrix[0]*fTrackInfo->dcosXVtx + rotMatrix[1]*fTrackInfo->dcosYVtx + rotMatrix[2]*fTrackInfo->dcosZVtx); dcosy_ideal = (rotMatrix[3]*fTrackInfo->dcosXVtx + rotMatrix[4]*fTrackInfo->dcosYVtx + rotMatrix[5]*fTrackInfo->dcosZVtx); dcosz_ideal = (rotMatrix[6]*fTrackInfo->dcosXVtx + rotMatrix[7]*fTrackInfo->dcosYVtx + rotMatrix[8]*fTrackInfo->dcosZVtx); fTrackInfo->azimuth = 180.*TMath::ATan2(-dcosx_ideal,dcosz_ideal)/TMath::Pi(); if(fTrackInfo->azimuth < 0.) fTrackInfo->azimuth += 360.; } if(ntpTrack->fit.ndof == 0 || ntpTrack->momentum.qp== 0.){ fGoodTrack = false; return; } //the value of trk.fidvtx.dz is not valid for events with the full detector //so hack up my own value of fiducialDz Float_t fullDetEndZ = 0.0594*485. + 1.; fTrackInfo->fiducialVtxDz = ntpTrack->fidvtx.dz; if(fEndPlane>248) fTrackInfo->fiducialVtxDz = TMath::Min(fullDetEndZ - fTrackInfo->vtxZ, fTrackInfo->vtxZ); fTrackInfo->fiducialVtxDr = ntpTrack->fidvtx.dr; fTrackInfo->fiducialEndDz = ntpTrack->fidend.dz; if(fEndPlane>248) fTrackInfo->fiducialEndDz = TMath::Min(fullDetEndZ - fTrackInfo->endZ, fTrackInfo->endZ); fTrackInfo->fiducialEndDr = ntpTrack->fidend.dr; fTrackInfo->uvAsymmetry = 1.; if(nVPlanes+nUPlanes>0) fTrackInfo->uvAsymmetry = TMath::Abs((1.*nUPlanes-1.*nVPlanes)/(1.*nVPlanes+1.*nUPlanes)); fTrackInfo->signalUseFraction = 0.; if(fEventInfo->pulseHeight>0.) fTrackInfo->signalUseFraction = fTrackInfo->pulseHeight/(fEventInfo->pulseHeight); fTrackInfo->planeUseFraction = 0.; if(fEventInfo->planes>0) fTrackInfo->planeUseFraction = 1.*fTrackInfo->trackLikePlanes/(1.*fEventInfo->planes); Int_t numStrips = 0; Int_t numDigits = 0; fTrackInfo->alternateChi2 = 0.; double deltaTPos = 0.; double stripPosError = 0.041/TMath::Sqrt(12.); fTrackInfo->impactParameter = 20.0; fTrackInfo->totalStrips = 0; //does the pathlength of the track make sense? if(fTrackInfo->length<50.){ //loop over all strips in the track numStrips = (int)ntpTrack->nstrip; numDoubleEndedStrips = 0; MSG("CondensedNtpModuleAtm", Msg::kDebug) << "start looping over track strips" << endl; for(Int_t ns = 0; ns < numStrips; ++ns){ //get the index for this strip index = ntpTrack->stp[ns]; //get the NtpSRStrip object ntpStrip = dynamic_cast<NtpSRStrip *>(stripArray->At(index)); plane = (int)ntpStrip->plane; strip = (int)ntpStrip->strip; numDigits = (int)ntpStrip->ndigit; xPos = ntpTrack->stpx[ns]; yPos = ntpTrack->stpy[ns]; if((int)ntpStrip->planeview == PlaneView::kU) deltaTPos = ntpTrack->stpu[ns]-ntpStrip->tpos; else if((int)ntpStrip->planeview == PlaneView::kV) deltaTPos = ntpTrack->stpv[ns]-ntpStrip->tpos; MSG("CondensedNtpModuleAtm", Msg::kDebug) << ntpTrack->stpu[ns] << " " << ntpStrip->tpos << " " << ntpTrack->stpv[ns] << " " << deltaTPos << endl; fTrackInfo->alternateChi2 += deltaTPos*deltaTPos/(stripPosError*stripPosError); if(numDigits >= 2){ ++numDoubleEndedStrips; radius = TMath::Sqrt(yPos*yPos + xPos*xPos); if(radius<fTrackInfo->impactParameter) fTrackInfo->impactParameter = TMath::Sqrt(yPos*yPos + xPos*xPos); //get values to determine the number of plances crossed by the track //in concentric circles in the detector from 0.3 m out to //4.0 m. the points come in order, so i dont have to test each one //to make sure it is before or after the previous points if(radius>0.3 && radius<1.0 && plane != prevPlane) ++fTrackInfo->planesIn10Meter; if(radius>0.3 && radius<1.5 && plane != prevPlane) ++fTrackInfo->planesIn15Meter; if(radius>0.3 && radius<2.0 && plane != prevPlane) ++fTrackInfo->planesIn20Meter; if(radius>0.3 && radius<2.5 && plane != prevPlane) ++fTrackInfo->planesIn25Meter; if(radius>0.3 && radius<3.0 && plane != prevPlane) ++fTrackInfo->planesIn30Meter; if(radius>0.3 && radius<3.5 && plane != prevPlane) ++fTrackInfo->planesIn35Meter; if(radius>0.3 && radius<4.0 && plane != prevPlane) ++fTrackInfo->planesIn40Meter; ++fTrackInfo->totalStrips; } // end if the strip had double sided readout prevPlane = plane; } // end loop over strips in the track fTrackInfo->twoEndStripFraction = 0.; if(numStrips>0){ fTrackInfo->twoEndStripFraction = (1.*numDoubleEndedStrips)/(1.*numStrips); fTrackInfo->alternateChi2 /= (1.*numStrips); } }//end if the track length is reasonable else fGoodTrack = false; MSG("CondensedNtpModuleAtm", Msg::kDebug) << "end fill track info - good track = " << (int)fGoodTrack << endl; return; }

void CondensedNtpModuleAtm::FillTrackMagneticBendingVariables (  NtpSRTrack *  ntpTrack, TClonesArray *  stripArray )  [private]

Definition at line 852 of file CondensedNtpModuleAtm.cxx. References NtpSRTrack::end, fGoodTrack, fTrackInfo, ANtpTrackInfoAtm::meanDistFromLinearFit, MSG, ANtpTrackInfoAtm::netDistFromLinearFit, NtpSRTrack::nstrip, ANtpTrackInfoAtm::rmsDistFromLinearFit, ANtpTrackInfoAtm::sagitta, NtpSRTrack::stp, NtpSRTrack::stpx, NtpSRTrack::stpy, NtpSRTrack::stpz, NtpSRTrack::vtx, NtpSRVertex::x, NtpSRVertex::y, NtpSRVertex::z, and ANtpTrackInfoAtm::zeroCurvatureChi2. Referenced by Ana(). { MSG("CondensedNtpModuleAtm", Msg::kDebug) << "in FillTrackMagneticBendingVariables" << endl; NtpSRStrip *ntpStrip = 0; //find the direction cosines for the straight line between the vertex and end points //of the track Float_t deltaX = ntpTrack->end.x - ntpTrack->vtx.x; Float_t deltaY = ntpTrack->end.y - ntpTrack->vtx.y; Float_t deltaZ = ntpTrack->end.z - ntpTrack->vtx.z; Float_t deltaS = TMath::Sqrt(deltaX*deltaX + deltaY*deltaY + deltaZ*deltaZ); Float_t dcosX = deltaX/deltaS; Float_t dcosY = deltaY/deltaS; Float_t dcosZ = deltaZ/deltaS; Float_t xPos = 0.; Float_t yPos = 0.; Float_t zPos = 0.; Float_t linXPos = 0.; Float_t linYPos = 0.; //loop over the strips in the track to find the sagitta fTrackInfo->sagitta = 0.; fTrackInfo->netDistFromLinearFit = 0.; fTrackInfo->meanDistFromLinearFit = 0.; fTrackInfo->rmsDistFromLinearFit = 0.; fTrackInfo->zeroCurvatureChi2 = 0.; if(deltaZ == 0. || deltaX == 0. || deltaY == 0.){ fGoodTrack = false; return; } TH1F dist("dist", "", 500, 0., 5.); Float_t curDistFromLinearFit = 0.; //loop over all strips in the track Int_t numStrips = (int)ntpTrack->nstrip; Int_t index = 0; double stripPosError = 0.041/TMath::Sqrt(12.); for(Int_t ns = 0; ns < numStrips; ++ns){ //get the index for this strip index = ntpTrack->stp[ns]; //get the NtpSRStrip object ntpStrip = dynamic_cast<NtpSRStrip *>(stripArray->At(index)); xPos = ntpTrack->stpx[ns]; yPos = ntpTrack->stpy[ns]; zPos = ntpTrack->stpz[ns]; linXPos = ntpTrack->vtx.x + (dcosX/dcosZ) * (zPos - ntpTrack->vtx.z); linYPos = ntpTrack->vtx.y + (dcosY/dcosZ) * (zPos - ntpTrack->vtx.z); deltaX = xPos - linXPos; deltaY = yPos - linYPos; fTrackInfo->zeroCurvatureChi2 += (deltaX*deltaX + deltaY*deltaY)/(stripPosError*stripPosError); curDistFromLinearFit = TMath::Sqrt(deltaX*deltaX + deltaY*deltaY); dist.Fill(curDistFromLinearFit); if(curDistFromLinearFit > fTrackInfo->sagitta) fTrackInfo->sagitta = curDistFromLinearFit; } if(numStrips>1) fTrackInfo->zeroCurvatureChi2 /= 1.*(numStrips-1); fTrackInfo->meanDistFromLinearFit = dist.GetMean(); fTrackInfo->rmsDistFromLinearFit = dist.GetRMS(); //now to figure out the net pt kick to the muon return; }

void CondensedNtpModuleAtm::FillTrackTimingVariables (  NtpSRTrack *  ntpTrack, TClonesArray *  stripArray )  [private]

Definition at line 710 of file CondensedNtpModuleAtm.cxx. References ANtpHeaderInfo::dataType, fGoodTrack, fHeaderInfo, FitMinimizingResidual(), ANtpTrackInfoAtm::flatLineChi2, fTrackInfo, GetTimeWeight(), ANtpTrackInfoAtm::invBeta, ANtpTrackInfoAtm::invBetaChi2, ANtpTrackInfo::length, MSG, NtpSRStrip::ndigit, NtpSRTrack::nstrip, NtpSRStrip::ph0, NtpSRStrip::ph1, NtpSRStrip::plane, ANtpHeaderInfo::run, NtpSRPulseHeight::sigcor, ANtpHeaderInfo::snarl, NtpSRTrack::stp, NtpSRTrack::stpds, NtpSRTrack::stpt0, NtpSRTrack::stpt1, NtpSRTrack::stpx, NtpSRTrack::stpy, NtpSRTrack::stpz, NtpSRStrip::strip, ANtpHeaderInfo::subRun, ANtpTrackInfoAtm::timeSlopeY, and ANtpTrackInfoAtm::timeSlopeYChi2. Referenced by Ana(). { MSG("CondensedNtpModuleAtm", Msg::kDebug) << "in FillTrackTimingVariables" << endl; NtpSRStrip *ntpStrip = 0; Double_t time[2] = {0.}; Double_t xPos = 0., yPos = 0., zPos = 0.; Double_t yPosition[5000] = {0.}; Double_t yPosition1[5000] = {0.}; Double_t weight[5000] = {0.}; Double_t weight2[5000] = {0.}; Double_t pathLength[5000] = {0.}; Double_t weight1[5000] = {0.}; Double_t timeAv[5000] = {0.}; Double_t timeAv1[5000] = {0.}; Int_t arrayCtr = 0, betaCtr = 0, actr = 0, plane = 0, strip = 0, index = 0; Double_t minTime = 1.e20; //loop over all strips in the track Int_t numStrips = (int)ntpTrack->nstrip; Int_t numDigits = 0; Double_t smear = 0.; MSG("CondensedNtpModuleAtm", Msg::kDebug) << "fill track timing" << endl; for(Int_t ns = 0; ns < numStrips; ++ns){ if(fHeaderInfo->dataType == (int)SimFlag::kMC) smear = gRandom->Gaus(0., 0.75); else smear = 0.; //get the index for this strip index = ntpTrack->stp[ns]; //get the NtpSRStrip object ntpStrip = dynamic_cast<NtpSRStrip *>(stripArray->At(index)); plane = (int)ntpStrip->plane; strip = (int)ntpStrip->strip; numDigits = (int)ntpStrip->ndigit; xPos = ntpTrack->stpx[ns]; yPos = ntpTrack->stpy[ns]; zPos = ntpTrack->stpz[ns]; if(numDigits == 2){ time[0] = ntpTrack->stpt0[ns]; time[1] = ntpTrack->stpt1[ns]; if(time[0] < minTime) minTime = time[0]; if(time[1] < minTime) minTime = time[1]; //have a point for each end of the strip timeAv[arrayCtr] = 1.e9*time[0] + smear; timeAv[arrayCtr+1] = 1.e9*time[1] + smear; pathLength[arrayCtr] = fTrackInfo->length - ntpTrack->stpds[ns]; pathLength[arrayCtr+1] = fTrackInfo->length - ntpTrack->stpds[ns]; //stuff to get timing information yPosition[arrayCtr] = yPos; weight[arrayCtr] = GetTimeWeight(ntpStrip->ph0.sigcor); if (weight[arrayCtr]>0.4) weight[arrayCtr]=0.4; yPosition[arrayCtr+1] = yPos; weight[arrayCtr+1] = GetTimeWeight(ntpStrip->ph1.sigcor); if (weight[arrayCtr+1]>0.4) weight[arrayCtr+1]=0.4; MSG("CondensedNtpModuleAtm", Msg::kDebug) << weight[arrayCtr] << " " << weight[arrayCtr+1] << endl; arrayCtr += 2; } // end if the strip had double sided readout } // end loop over strips in the track //subtract the minimum time and copy the weights into another array to do the fit to //a flat line later for(Int_t ti = 0; ti< arrayCtr; ++ti){ timeAv[ti] -= 1.e9*minTime; weight2[ti] = weight[ti]; } //cout << "got strip info" << endl; //array for fit parameters Double_t param[] = {-10000., -10000.}; Double_t eparam[] = {-10000., -10000.}; //find the inverse beta value betaCtr = arrayCtr; FitMinimizingResidual(betaCtr, pathLength, timeAv, weight, param, eparam); fTrackInfo->invBeta = param[1]*0.3; if(betaCtr > 2) fTrackInfo->invBetaChi2 = eparam[0] /(1.*(betaCtr-2)); //find the chi^2 for a fit to a straight line so you //can see if the sloped fit is really better FitMinimizingResidual(betaCtr, pathLength, timeAv, weight2, param, eparam, false); if(betaCtr > 1) fTrackInfo->flatLineChi2 = eparam[0] / (1.*betaCtr-1); //now loop over the points in the array and fill new arrays //containing only those points that have a nonzero weight actr = 0; MSG("CondensedNtpModuleAtm", Msg::kDebug) << fHeaderInfo->run << " " << fHeaderInfo->subRun << " " << fHeaderInfo->snarl << " " << arrayCtr << endl; for(Int_t ii = 0; ii < arrayCtr; ++ii){ MSG("CondensedNtpModuleAtm", Msg::kDebug) << " " << ii << "/" << arrayCtr << " " << yPosition[ii] << " " << weight[ii] << " " << timeAv[ii] << " " << actr << endl; if(weight[ii]>0.){ yPosition1[actr] = yPosition[ii]; weight1[actr] = weight[ii]; timeAv1[actr] = timeAv[ii]; ++actr; } }//end loop to select good points along the track //get the fit for the time in the y direction FitMinimizingResidual(actr, yPosition1, timeAv1, weight1, param, eparam); fTrackInfo->timeSlopeY = param[1]; fTrackInfo->timeSlopeYChi2 = eparam[0]; if(actr > 2) fTrackInfo->timeSlopeYChi2 /= 1.*(actr-2); MSG("CondensedNtpModuleAtm", Msg::kDebug) << "actr = " << actr << "/" << arrayCtr << " goodtrack = " << (int)fGoodTrack << endl; if(1.*actr < 0.5*arrayCtr) fGoodTrack = false; return; }

void CondensedNtpModuleAtm::FindChi2 (  Double_t *  x, Double_t *  y, Double_t *  weight, Int_t  nPoints, Double_t &  a1, Double_t &  a2, Double_t &  chiSq )  [private]

Definition at line 1277 of file CondensedNtpModuleAtm.cxx. { double vtxX = x[0]; double vtxY = y[0]; double endX = x[nPoints-1]; double endY = y[nPoints-1]; a2 = endY-vtxY; if(TMath::Abs(vtxX-endX)>0.) a2 /= (endX-vtxX); else{ cout << "cannot find slope for n = " << nPoints << " points " << endl; for(Int_t i = 0; i<nPoints; ++i) cout << x[i] << " " << y[i] << " " << weight[i] << endl; } a1 = endY-a2*endX; //find the chi^2 value for the fit chiSq = 0.; Double_t chi = 0.; for(Int_t j = 0; j<nPoints; j++){ chi = (y[j] - (a1 + a2*x[j]))/weight[j]; // cout << y[j] << " " << a1+a2*x[j] << " " << chiSq << endl; chiSq += chi*chi; } // cout << "chi^2 = " << chiSq << " " << a1 << " " << a2 << endl; //get the chiSq per degree of freedom if(nPoints>2) chiSq /= 1.*(nPoints-2); return; }

void CondensedNtpModuleAtm::FindLinearFit (  Double_t *  x, Double_t *  y, Double_t *  weight, Int_t  nPoints, Double_t &  a1, Double_t &  a2, Double_t &  chiSq )  [private]

Definition at line 1179 of file CondensedNtpModuleAtm.cxx. { //use method of determinants to do the least squares fit //this is from bevington chapter 6 // //f_1(x) = 1 //f_2(x) = x //need to find the sums //SIGMA y_i*f_1(x_i)/weight_i //SIGMA y_i*f_2(x_i)/weight_i //SIGMA f_1(x_i)*f_1(x_i)/weight_i //SIGMA f_1(x_i)*f_2(x_i)/weight_i //SIGMA f_2(x_i)*f_1(x_i)/weight_i //SIGMA f_2(x_i)*f_2(x_i)/weight_i // //some of these are the same, ie f_2*f_1 = f_1*f_2 so only need //3 instead of 4 variables Double_t yf1divWeight = 0.; Double_t yf2divWeight = 0.; Double_t f1f1divWeight = 0.; Double_t f1f2divWeight = 0.; Double_t f2f2divWeight = 0.; for(Int_t i=0; i<nPoints; i++){ // cout << y[i] << " " << x[i] << " " << weight[i] << endl; yf1divWeight += y[i]/(weight[i]*weight[i]); yf2divWeight += (y[i]*x[i])/(weight[i]*weight[i]); f1f1divWeight += 1./(weight[i]*weight[i]); f1f2divWeight += x[i]/(weight[i]*weight[i]); f2f2divWeight += (x[i]*x[i])/(weight[i]*weight[i]); } //make a bunch of matricies to hold these values TMatrix deltaMatrix(2,2); TMatrix a1Matrix(2,2); TMatrix a2Matrix(2,2); //fill the matricies deltaMatrix(0,0) = f1f1divWeight; deltaMatrix(0,1) = f1f2divWeight; deltaMatrix(1,0) = f1f2divWeight; deltaMatrix(1,1) = f2f2divWeight; a1Matrix(0,0) = yf1divWeight; a1Matrix(0,1) = f1f2divWeight; a1Matrix(1,0) = yf2divWeight; a1Matrix(1,1) = f2f2divWeight; a2Matrix(0,0) = f1f1divWeight; a2Matrix(0,1) = yf1divWeight; a2Matrix(1,0) = f1f2divWeight; a2Matrix(1,1) = yf2divWeight; a1 = -10000.; a2 = -10000.; Double_t delta = deltaMatrix.Determinant(); // cout << "delta = " << delta << endl; if(delta != 0.){ a1 = a1Matrix.Determinant()/delta; a2 = a2Matrix.Determinant()/delta; } //find the chi^2 value for the fit chiSq = 0.; Double_t chi = 0.; for(Int_t j = 0; j<nPoints; j++){ chi = (y[j] - (a1 + a2*x[j]))/weight[j]; // cout << chi*weight[j] << endl; chiSq += chi*chi; } // cout << "chi^2 = " << chiSq << " " << a1 << " " << a2 << endl; //get the chiSq per degree of freedom chiSq /= (1.*nPoints - 2.); return; }

void CondensedNtpModuleAtm::FindMeanAndRMS (  double *  x, double *  weight, int  nPoints, double &  mean, double &  rms )  [private]

Definition at line 1366 of file CondensedNtpModuleAtm.cxx. { // cout << "rms" << endl; double sumXXW = 0.; double sumXW = 0.; double sumW = 0.; for(Int_t i = 0; i<nPoints; ++i){ sumXXW += x[i]*x[i]*weight[i]; sumXW += x[i]*weight[i]; sumW += weight[i]; // cout << sumXW << " " << sumXXW << " " << sumW << " " // << x[i] << " " << weight[i] << endl; } if(sumW > 0. && nPoints>1){ mean = sumXW/sumW; if(sumXXW >=0.) rms = TMath::Sqrt((sumXXW/sumW - mean*mean) *(nPoints*1./(1.*nPoints-1.))); } else{ rms = 0.; mean = x[0]; } return; }

void CondensedNtpModuleAtm::FindStraightLineDeviation (  Double_t *  x, Double_t *  y, Int_t  nPoints, Double_t &  a1, Double_t &  a2, Double_t &  deviation )  [private]

Definition at line 1325 of file CondensedNtpModuleAtm.cxx. { double vtxX = x[0]; double vtxY = y[0]; double endX = x[nPoints-1]; double endY = y[nPoints-1]; a2 = endY-vtxY; if(TMath::Abs(vtxX-endX)>0.) a2 /= (endX-vtxX); else{ cout << "cannot find slope for n = " << nPoints << " points " << endl; for(Int_t i = 0; i<nPoints; ++i) cout << x[i] << " " << y[i] << endl; } a1 = endY-a2*endX; //find the chi^2 value for the fit deviation = 0.; Double_t dev = 0.; for(Int_t j = 0; j<nPoints; j++){ dev = (y[j] - (a1 + a2*x[j])); // cout << y[j] << " " << a1+a2*x[j] << " " << chiSq << endl; deviation += dev*dev; } return; }

void CondensedNtpModuleAtm::FitMinimizingResidual (  Int_t &  n, Double_t *  x, Double_t *  y, Double_t *  w, Double_t *  param, Double_t *  eparam, bool  findSlope = true )  [private]

Definition at line 1118 of file CondensedNtpModuleAtm.cxx. References LinearFit_Weighted(), and WeightedAverage(). Referenced by FillTrackTimingVariables(). { // cout << "in LinearFit_MinimizeResid" << endl; Int_t nremove = 0; Double_t maxresid = 11.; Double_t resid = 0.; Double_t parm[2] = {0.}; Double_t eparm[2] = {0.}; Double_t timefitchi2 = 0.; while(maxresid>10. &&n-nremove>4 &&(1.*nremove)/(1.*n)<0.2){ timefitchi2=0.; if(findSlope) LinearFit_Weighted(n,x,y,w,parm,eparm); else WeightedAverage(n,y,w,parm,eparm); maxresid = 0.; Int_t imaxresid=0; for (Int_t i=0; i<n; i++) { if (w[i]>0.) { if(findSlope) resid = (parm[1]*x[i]+parm[0]-y[i]); else resid = parm[0] - y[i]; timefitchi2 += w[i]*resid*resid; if (TMath::Abs(resid)>maxresid) { maxresid = TMath::Abs(resid); imaxresid = i; } } } if (maxresid>10.) { w[imaxresid] = 0.; nremove++; } }//end loop to minimize residuals //get the chi^2 timefitchi2 = 0.; for (Int_t i=0; i<n; i++) { if (w[i]>0.) { if(findSlope) resid = (parm[1]*x[i]+parm[0]-y[i]); else resid = parm[0] - y[i]; timefitchi2 += w[i]*resid*resid; } } n -= nremove; param[0] = parm[0]; param[1] = parm[1]; eparam[0] = timefitchi2; eparam[1] = eparm[1]; return; }

Double_t CondensedNtpModuleAtm::GetTimeWeight (  Float_t  ph  )  [private]

Definition at line 1000 of file CondensedNtpModuleAtm.cxx. References ArrivalTime_Weight(). Referenced by FillTrackTimingVariables(). { // cout << "in GetTimeWeight" << endl; Double_t npe = ph/60.; if (npe<1.) npe=1.; // cout << ArrivalTime_Weight(npe) << endl; return ArrivalTime_Weight(npe); }

void CondensedNtpModuleAtm::Help (   )  [virtual]

Implement to spew some useful help to cout Reimplemented from JobCModule. Definition at line 410 of file CondensedNtpModuleAtm.cxx. References MSG. { MSG("JobC", Msg::kInfo) << "NearElectronicsCheck::Help\n" <<"CondensedNtpModuleAtm is a module which demultiplexes events " <<"in the far detector." << endl; }

void CondensedNtpModuleAtm::LinearFit_Weighted (  const Int_t  n, const Double_t *  x, const Double_t *  y, const Double_t *  w, Double_t *  parm, Double_t *  eparm )  [private]

Definition at line 1050 of file CondensedNtpModuleAtm.cxx. Referenced by FitMinimizingResidual(). { // cout << "in LinearFit_Weighted with errors" << endl; Double_t sumx=0.; Double_t sumx2=0.; Double_t sumy=0.; Double_t sumy2=0.; Double_t sumxy=0.; Double_t sumw=0.; for(Int_t i=0; i<n; ++i){ sumx += x[i]*w[i]; sumx2 += x[i]*x[i]*w[i]; sumy += y[i]*w[i]; sumy2 += y[i]*y[i]*w[i]; sumxy += x[i]*y[i]*w[i]; sumw += w[i]; // cout << x[i] << " " << y[i] << " " << w[i] << endl; } // cout << sumx2 << " " << sumw << " " << sumx << " " << sumxy << " " << sumy // << endl; if(sumx2*sumw-sumx*sumx>0.&&sumx2<1.e5){ parm[0] = (sumy*sumx2-sumx*sumxy)/(sumx2*sumw-sumx*sumx); parm[1] = (sumxy*sumw-sumx*sumy)/(sumx2*sumw-sumx*sumx); eparm[0] = TMath::Sqrt(sumx2*(sumx2*sumw-sumx*sumx))/(sumx2*sumw-sumx*sumx); eparm[1] = TMath::Sqrt(sumx2*sumw-sumx*sumx)/(sumx2*sumw-sumx*sumx); } return; }

void CondensedNtpModuleAtm::ResetTreeVariables (   )  [private]

Definition at line 989 of file CondensedNtpModuleAtm.cxx. References fEventInfo, fGoodTrack, fTrackInfo, fTruthInfo, ANtpTruthInfoAtm::Reset(), ANtpTrackInfoAtm::Reset(), and ANtpEventInfo::Reset(). Referenced by Ana(). { fEventInfo->Reset(); fTrackInfo->Reset(); fTruthInfo->Reset(); fGoodTrack = false; return; }

void CondensedNtpModuleAtm::WeightedAverage (  const Int_t  n, const Double_t *  x, const Double_t *  w, Double_t *  parm, Double_t *  eparm )  [private]

Definition at line 1091 of file CondensedNtpModuleAtm.cxx. Referenced by FitMinimizingResidual(). { // cout << "in LinearFit_Weighted with errors" << endl; double sumWeights = 0.; double sumXWeights = 0.; //find the weighted average of the points for(int i = 0; i < n; ++i){ sumXWeights += x[i]*w[i]; sumWeights += w[i]; } if(sumWeights > 0.) parm[0] = sumXWeights/sumWeights; else parm[0] = -9999.; parm[1] = 0.; eparm[0] = 0.; eparm[1] = 0.; return; }

---

Member Data Documentation

Double_t CondensedNtpModuleAtm::fBeamHPos [private]

Definition at line 72 of file CondensedNtpModuleAtm.h.

ANtpBeamInfo* CondensedNtpModuleAtm::fBeamInfo [private]

Definition at line 88 of file CondensedNtpModuleAtm.h. Referenced by Ana(), and BeginJob().

std::string CondensedNtpModuleAtm::fBeamPath [private]

Definition at line 61 of file CondensedNtpModuleAtm.h. Referenced by BeginJob(), and Config().

Double_t CondensedNtpModuleAtm::fBeamTimeSec [private]

Definition at line 75 of file CondensedNtpModuleAtm.h.

TTree* CondensedNtpModuleAtm::fBeamTree [private]

Definition at line 65 of file CondensedNtpModuleAtm.h.

std::string CondensedNtpModuleAtm::fBeamTreeName [private]

Definition at line 60 of file CondensedNtpModuleAtm.h. Referenced by BeginJob(), and Config().

Double_t CondensedNtpModuleAtm::fBeamVPos [private]

Definition at line 73 of file CondensedNtpModuleAtm.h.

Double_t CondensedNtpModuleAtm::fColDenCosZen[1500] [private]

Definition at line 82 of file CondensedNtpModuleAtm.h. Referenced by BeginJob().

Int_t CondensedNtpModuleAtm::fCtr [private]

Definition at line 94 of file CondensedNtpModuleAtm.h.

Double_t CondensedNtpModuleAtm::fCurrentPOT [private]

Definition at line 71 of file CondensedNtpModuleAtm.h.

Int_t CondensedNtpModuleAtm::fDataType [private]

Definition at line 68 of file CondensedNtpModuleAtm.h. Referenced by BeginJob(), and Config().

Double_t CondensedNtpModuleAtm::fDensity[1500] [private]

Definition at line 83 of file CondensedNtpModuleAtm.h. Referenced by BeginJob(), and FillTrackInformation().

Int_t CondensedNtpModuleAtm::fDetector [private]

Definition at line 66 of file CondensedNtpModuleAtm.h. Referenced by Ana(), Config(), and FillTrackInformation().

Int_t CondensedNtpModuleAtm::fEndPlane [private]

Definition at line 69 of file CondensedNtpModuleAtm.h. Referenced by Ana(), Config(), and FillTrackInformation().

ANtpEventInfo* CondensedNtpModuleAtm::fEventInfo [private]

Definition at line 86 of file CondensedNtpModuleAtm.h. Referenced by Ana(), BeginJob(), FillEventInformation(), FillTrackInformation(), and ResetTreeVariables().

Int_t CondensedNtpModuleAtm::fFailDeMux [private]

Definition at line 78 of file CondensedNtpModuleAtm.h. Referenced by Ana(), and BeginJob().

TTree* CondensedNtpModuleAtm::fFailTree [private]

Definition at line 64 of file CondensedNtpModuleAtm.h. Referenced by Ana(), BeginJob(), and EndJob().

std::string CondensedNtpModuleAtm::fFileName [private]

Definition at line 57 of file CondensedNtpModuleAtm.h. Referenced by BeginJob(), and Config().

Int_t CondensedNtpModuleAtm::fGoodTrack [private]

Definition at line 67 of file CondensedNtpModuleAtm.h. Referenced by Ana(), FillTrackInformation(), FillTrackMagneticBendingVariables(), FillTrackTimingVariables(), and ResetTreeVariables().

ANtpHeaderInfo* CondensedNtpModuleAtm::fHeaderInfo [private]

Definition at line 87 of file CondensedNtpModuleAtm.h. Referenced by Ana(), BeginJob(), and FillTrackTimingVariables().

Double_t CondensedNtpModuleAtm::fHornCurrent [private]

Definition at line 70 of file CondensedNtpModuleAtm.h.

ANtpInfoObjectFillerBeam* CondensedNtpModuleAtm::fInfoFiller [private]

Definition at line 92 of file CondensedNtpModuleAtm.h. Referenced by Ana(), FillEventInformation(), FillMCInformation(), and FillTrackInformation().

Int_t CondensedNtpModuleAtm::fIsMultiple [private]

Definition at line 76 of file CondensedNtpModuleAtm.h. Referenced by Ana(), and BeginJob().

Int_t CondensedNtpModuleAtm::fMajorRelease [private]

Definition at line 79 of file CondensedNtpModuleAtm.h. Referenced by Ana(), and Config().

bool CondensedNtpModuleAtm::fNewMC [private]

Definition at line 84 of file CondensedNtpModuleAtm.h. Referenced by Config().

TFile* CondensedNtpModuleAtm::fNtpFile [private]

Definition at line 62 of file CondensedNtpModuleAtm.h. Referenced by BeginJob(), and EndJob().

TTree* CondensedNtpModuleAtm::fNtuple [private]

Definition at line 63 of file CondensedNtpModuleAtm.h. Referenced by Ana(), BeginJob(), and EndJob().

Double_t CondensedNtpModuleAtm::fPathAzimuthDeg[1500] [private]

Definition at line 81 of file CondensedNtpModuleAtm.h. Referenced by BeginJob().

Double_t CondensedNtpModuleAtm::fPathCosZenDeg[1500] [private]

Definition at line 80 of file CondensedNtpModuleAtm.h. Referenced by BeginJob().

std::string CondensedNtpModuleAtm::fRockMapFileName [private]

Definition at line 58 of file CondensedNtpModuleAtm.h. Referenced by BeginJob(), and Config().

Double_t CondensedNtpModuleAtm::fTargetPos [private]

Definition at line 74 of file CondensedNtpModuleAtm.h.

ANtpTrackInfoAtm* CondensedNtpModuleAtm::fTrackInfo [private]

Definition at line 89 of file CondensedNtpModuleAtm.h. Referenced by Ana(), BeginJob(), FillTrackInformation(), FillTrackMagneticBendingVariables(), FillTrackTimingVariables(), and ResetTreeVariables().

std::string CondensedNtpModuleAtm::fTreeName [private]

Definition at line 59 of file CondensedNtpModuleAtm.h. Referenced by BeginJob(), and Config().

ANtpTruthInfoAtm* CondensedNtpModuleAtm::fTruthInfo [private]

Definition at line 90 of file CondensedNtpModuleAtm.h. Referenced by BeginJob(), FillMCInformation(), and ResetTreeVariables().

Int_t CondensedNtpModuleAtm::fValidPlaneFail [private]

Definition at line 77 of file CondensedNtpModuleAtm.h. Referenced by Ana(), and BeginJob().

---

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

• CondensedNtpModuleAtm.h
• CondensedNtpModuleAtm.cxx

---