CondensedNtpModuleAtm Class Reference

#include <CondensedNtpModuleAtm.h>

Inheritance diagram for CondensedNtpModuleAtm:

JobCModule List of all members.

Public Member Functions

 CondensedNtpModuleAtm ()
virtual ~CondensedNtpModuleAtm ()
void BeginJob ()
JobCResult Ana (const MomNavigator *mom)
const RegistryDefaultConfig () 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
ANtpEventInfofEventInfo
ANtpHeaderInfofHeaderInfo
ANtpBeamInfofBeamInfo
ANtpTrackInfoAtmfTrackInfo
ANtpTruthInfoAtmfTruthInfo
ANtpInfoObjectFillerBeamfInfoFiller
Int_t fCtr

Detailed Description

Definition at line 41 of file CondensedNtpModuleAtm.h.


Constructor & Destructor Documentation

CondensedNtpModuleAtm::CondensedNtpModuleAtm (  ) 

Definition at line 70 of file CondensedNtpModuleAtm.cxx.

References Msg::kDebug, and MSG.

00070                                              : 
00071     fFileName("analysisNtuple.root"),
00072     fTreeName("analysisNtuple"),
00073     fNtpFile(0),
00074     fNtuple(0),
00075     fBeamTree(0),
00076     fDetector(1),
00077     fEndPlane(485),
00078     fEventInfo(0),
00079     fHeaderInfo(0),
00080     fTrackInfo(0),
00081     fTruthInfo(0),
00082     fCtr(0)
00083 {
00084     MSG("JobC", Msg::kDebug) << "CondensedNtpModuleAtm::Constructor" << endl;
00085 
00086     fHeaderInfo = new ANtpHeaderInfo();
00087     fBeamInfo = new ANtpBeamInfo();
00088     fEventInfo = new ANtpEventInfo();
00089     fTrackInfo = new ANtpTrackInfoAtm();
00090     fTruthInfo = new ANtpTruthInfoAtm();
00091     
00092     fInfoFiller = new ANtpInfoObjectFillerBeam();
00093 
00094 }

CondensedNtpModuleAtm::~CondensedNtpModuleAtm (  )  [virtual]

Definition at line 97 of file CondensedNtpModuleAtm.cxx.

References fBeamInfo, fEventInfo, fHeaderInfo, fTrackInfo, fTruthInfo, Msg::kDebug, and MSG.

00098 {
00099 
00100     MSG("JobC", Msg::kDebug) << "CondensedNtpModuleAtm::Destructor" << endl;
00101 
00102     if(fHeaderInfo) delete fHeaderInfo; 
00103     if(fBeamInfo) delete fBeamInfo;     
00104     if(fEventInfo) delete fEventInfo;   
00105     if(fTrackInfo) delete fTrackInfo;   
00106     if(fTruthInfo) delete fTruthInfo;   
00107     
00108 }


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 det, ANtpHeaderInfo::detector, fBeamInfo, fCtr, fDataType, fDetector, fEndPlane, fEventInfo, fFailDeMux, fFailTree, fGoodTrack, fHeaderInfo, ANtpInfoObjectFillerBeam::FillBeamInformation(), FillEventInformation(), ANtpInfoObjectFiller::FillHeaderInformation(), FillMCInformation(), FillTrackInformation(), FillTrackMagneticBendingVariables(), FillTrackTimingVariables(), fInfoFiller, fIsMultiple, fMajorRelease, fNewMC, fNtuple, fTrackInfo, fValidPlaneFail, ANtpRecoNtpManipulator::GetDmxStatusInfo(), ANtpEventManipulator::GetEvent(), ANtpRecoNtpManipulator::GetEventManipulator(), MomNavigator::GetFragment(), ANtpRecoNtpManipulator::GetSnarlEventSummary(), ANtpRecoNtpManipulator::GetStripArray(), VldContext::GetTimeStamp(), RecRecordImp< T >::GetVldContext(), ANtpEventInfo::index, NtpSRDmxStatus::ismultimuon, SimFlag::kData, Msg::kDebug, Detector::kFar, SimFlag::kMC, Detector::kNear, JobCResult::kPassed, Msg::kWarning, MSG, NtpSREventSummary::nevent, NtpSRDmxStatus::nonphysicalfail, NtpSREventSummary::ntrack, ANtpTrackInfo::planes, 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.

00191 {
00192   MSG("CondensedNtpModuleAtm", Msg::kDebug) << "start ana" << endl;
00193   
00194   JobCResult result(JobCResult::kPassed);
00195   
00196   //get the records from MOM
00197   assert(mom);
00198     
00199   NtpSRRecord *record = dynamic_cast<NtpSRRecord *>(mom->GetFragment("NtpSRRecord"));
00200   NtpStRecord *stRecord = dynamic_cast<NtpStRecord *>(mom->GetFragment("NtpStRecord"));
00201   
00202   if(!record && !stRecord){
00203     MSG("CondensedNtpModuleAtm", Msg::kWarning) << "Could not get NtpSR or NtpSt Record from MOM" << endl;
00204     result.SetWarning().SetFailed();
00205     return result;
00206   }//end if no NtpSRRecord
00207   
00208   NtpMCRecord *mcRecord = 0;
00209   NtpTHRecord *thRecord = 0;
00210   
00211   mcRecord = dynamic_cast<NtpMCRecord *>(mom->GetFragment("NtpMCRecord"));
00212   thRecord = dynamic_cast<NtpTHRecord *>(mom->GetFragment("NtpTHRecord"));              
00213   
00214   //get the record for the beam tree
00215   NtpBDLiteRecord *bdRecord = dynamic_cast<NtpBDLiteRecord *>(mom->GetFragment("NtpBDLiteRecord"));
00216   if(!bdRecord) MSG("CondensedNtpModuleAtm", Msg::kDebug) << "NO BEAM RECORD" << endl;
00217 
00218   //fill the header information
00219   Detector::Detector_t det = Detector::kFar;
00220   if(fDetector<1) det = Detector::kNear;
00221   //get a sim flag
00222   SimFlag::SimFlag_t dataType = SimFlag::kData;
00223   if(fDataType) dataType = SimFlag::kMC;
00224 
00225   //instantiate a NtpHelper object to help you get the info you want
00226   ANtpRecoNtpManipulator *ntpManipulator = 0;
00227   if(record) 
00228     ntpManipulator = new ANtpRecoNtpManipulator(record, mcRecord, thRecord);
00229   else if(stRecord) 
00230     ntpManipulator = new ANtpRecoNtpManipulator(stRecord);
00231     
00232   VldTimeStamp timeStamp = stRecord ? stRecord->GetVldContext()->GetTimeStamp() : record->GetVldContext()->GetTimeStamp();
00233 
00234   MSG("CondensedNtpModuleAtm", Msg::kDebug) << timeStamp << endl;
00235 
00236 //   TStopwatch stopwatch;
00237 //   stopwatch.Start();
00238 
00239   fInfoFiller->SetStripArray(ntpManipulator->GetStripArray());
00240   fInfoFiller->FillHeaderInformation(ntpManipulator, det, dataType, fHeaderInfo);
00241 
00242 //   stopwatch.Stop();
00243 //   MSG("CondensedNtpModuleAtm", Msg::kDebug) << "took " << stopwatch.RealTime()
00244 //                                         << " to fill header" << endl;
00245 //   stopwatch.Reset();
00246 
00247 //   stopwatch.Start();
00248 
00249   if(fHeaderInfo->year > 2004){
00250     if(bdRecord) fInfoFiller->FillBeamInformation(bdRecord, fBeamInfo);
00251     else fInfoFiller->FillBeamInformation(timeStamp, det, dataType, fBeamInfo);
00252   }
00253 
00254 //   stopwatch.Stop();
00255 //   MSG("CondensedNtpModuleAtm", Msg::kDebug) << "took " << stopwatch.RealTime()
00256 //                                         << " to fill beam" << endl;
00257 //   stopwatch.Reset();
00258 
00259   //set up which flags you want to use to determine the primary shower or track
00260   //a value of 0 for a flag means it will not be used
00261   ntpManipulator->SetPrimaryTrackCriteria(0,1,0); // nplanes, length, total pulse height
00262   ntpManipulator->SetPrimaryShowerCriteria(0,1); // nplanes, total pulse height
00263   
00264   //check that the current event passed demuxing
00265   bool passedDeMux = true; 
00266   fGoodTrack = true;
00267   fIsMultiple = 0;
00268   fValidPlaneFail = 0;
00269   fFailDeMux = 0;
00270 
00271   MSG("CondensedNtpModuleAtm", Msg::kDebug) << "run,snarl = " << fHeaderInfo->run 
00272                                             << "," << fHeaderInfo->snarl << " " 
00273                                             << ntpManipulator->GetSnarlEventSummary().ntrack 
00274                                             << " " << passedDeMux << endl;
00275 
00276   if(fHeaderInfo->detector == Detector::kFar){
00277 
00278     if(ntpManipulator->GetDmxStatusInfo().ismultimuon){
00279       passedDeMux = false;
00280       fIsMultiple = 1;
00281     }
00282     if(ntpManipulator->GetDmxStatusInfo().validplanesfail){
00283       passedDeMux = false;
00284       fValidPlaneFail = 1;
00285     }
00286     if(fMajorRelease < 18){
00287 
00288       if(ntpManipulator->GetDmxStatusInfo().vertexplanefail
00289          || ntpManipulator->GetDmxStatusInfo().nonphysicalfail) passedDeMux = false;
00290       
00291       UShort_t uValidPlanes = ntpManipulator->GetDmxStatusInfo().uvalidplanes;
00292       UShort_t uStrayPlanes = ntpManipulator->GetDmxStatusInfo().ustrayplanes;
00293       UShort_t vValidPlanes = ntpManipulator->GetDmxStatusInfo().vvalidplanes;
00294       UShort_t vStrayPlanes = ntpManipulator->GetDmxStatusInfo().vstrayplanes;
00295       Float_t ufrac = 0., vfrac = 0.;
00296     
00297       if(uValidPlanes>0) 
00298         ufrac = (1.*uStrayPlanes)/(1.*uValidPlanes);
00299       else ufrac = 1.;
00300       if(vValidPlanes>0) 
00301         vfrac = (1.*vStrayPlanes)/(1.*vValidPlanes);
00302       else vfrac = 1.;
00303       if(ufrac>=0.45) passedDeMux = false;
00304       else if( vfrac>=0.45 ) passedDeMux = false;
00305       else if(uStrayPlanes>=3&&uValidPlanes<20) passedDeMux = false;
00306       else if(vStrayPlanes>=3&&vValidPlanes<20) passedDeMux = false;
00307       else if(uStrayPlanes>=4&&uValidPlanes>=20 
00308               && uValidPlanes<=40) passedDeMux = false;
00309       else if(vStrayPlanes>=4&&vValidPlanes>=20 
00310               && vValidPlanes<=40) passedDeMux = false;
00311       else if(ufrac>=0.1 && uValidPlanes>=40) passedDeMux = false;
00312       else if(vfrac>=0.1 && vValidPlanes>=40) passedDeMux = false;
00313     
00314     }//end if before switch over to cambridge demuxing
00315   }//end if far detector
00316 
00317   //if the event isnt demuxed well, then forget it also forget it
00318   //if there arent any tracks in the event
00319   if(!passedDeMux || ntpManipulator->GetSnarlEventSummary().ntrack==0 || fEndPlane<248){
00320     fFailTree->Fill();
00321     return result;
00322   }
00323 
00324   //set up some pointers to typical NtpSR objects
00325   NtpSRTrack *ntpTrack = 0;
00326   NtpSREvent *ntpEvent = 0;
00327   NtpMCTruth *ntpMCTruth = 0;
00328   NtpMCStdHep *ntpMCStdHep = 0;
00329   NtpTHEvent *ntpTHEvent = 0;
00330   
00331   //loop over the events in the snarl
00332   for(Int_t i = 0; i < ntpManipulator->GetSnarlEventSummary().nevent; ++i){
00333     
00334     MSG("CondensedNtpModuleAtm", Msg::kDebug) << "on event " << i << endl;
00335     
00336     ResetTreeVariables();
00337     
00338     fEventInfo->index = i;
00339 
00340 //     stopwatch.Start();
00341     //get event i.  the call to ANtpRecoNtpManipulator::GetEventInSnarl sets the 
00342     //NtpSREvent data member in that object to the current event so that
00343     //you can later call for the primary track, shower, etc.
00344     ntpManipulator->GetEventManipulator()->SetEventInSnarl(i);          
00345     ntpEvent = ntpManipulator->GetEventManipulator()->GetEvent();               
00346     if(ntpEvent) FillEventInformation(ntpManipulator,ntpEvent);
00347 //     stopwatch.Stop();
00348 //     MSG("CondensedNtpModuleAtm", Msg::kInfo) << "took " << stopwatch.RealTime()
00349 //                                           << " to fill event" << endl;
00350 //     stopwatch.Reset();
00351     
00352 //     stopwatch.Start();
00353     //get the primary track for the event - if no track is present it 
00354     //returns 0
00355     fGoodTrack = false;
00356     ntpTrack = ntpManipulator->GetEventManipulator()->GetPrimaryTrack();                
00357     if(ntpTrack){
00358       fGoodTrack = true;
00359       FillTrackInformation(ntpTrack, ntpManipulator->GetStripArray(), 
00360                            ntpManipulator->GetCosmicRayInfo().azimuth,
00361                            ntpManipulator->GetCosmicRayInfo().ra,
00362                            ntpManipulator->GetCosmicRayInfo().dec);
00363       FillTrackMagneticBendingVariables(ntpTrack, ntpManipulator->GetStripArray());
00364       if(fDetector>0&&fGoodTrack) FillTrackTimingVariables(ntpTrack, ntpManipulator->GetStripArray());
00365     } 
00366 //     stopwatch.Stop();
00367 //     MSG("CondensedNtpModuleAtm", Msg::kInfo) << "took " << stopwatch.RealTime()
00368 //                                           << " to fill track" << endl;
00369 //     stopwatch.Reset();
00370 
00371 //     stopwatch.Start();
00372     // Get best neu match from Truthhelper
00373     ntpTHEvent = ntpManipulator->GetMCManipulator()->GetNtpTHEvent(i);
00374     MSG("CondensedNtpModuleAtm", Msg::kDebug) << "start filling mc info " << ntpTHEvent->neumc << endl;
00375     if(ntpTHEvent) ntpMCTruth = ntpManipulator->GetMCManipulator()->GetNtpMCTruth(ntpTHEvent->neumc);
00376     if(ntpMCTruth){
00377       //get the NtpMCStdHep object - use the fact that i put the values i want in the parent 
00378       //place, which means i want the 0th entry
00379       if(fNewMC) ntpMCStdHep = ntpManipulator->GetMCManipulator()->GetNtpMCStdHep(0);
00380       FillMCInformation(ntpMCTruth, ntpMCStdHep);       
00381     }
00382 //     stopwatch.Stop();
00383 //     MSG("CondensedNtpModuleAtm", Msg::kInfo) << "took " << stopwatch.RealTime()
00384 //                                           << " to fill truth" << endl;
00385 //     stopwatch.Reset();
00386 
00387     if(fEventInfo->totalStrips<fTrackInfo->totalStrips) fGoodTrack = false;
00388     if(fTrackInfo->planes<1) fGoodTrack = false;
00389     if(fTrackInfo->totalStrips<1) fGoodTrack = false;
00390     
00391     MSG("CondensedNtpModuleAtm", Msg::kDebug) << "strips = " << fEventInfo->totalStrips 
00392                                               << "/" << fTrackInfo->totalStrips << " planes = " 
00393                                               << fTrackInfo->planes << " good track = " 
00394                                               << (int)fGoodTrack << endl;
00395     
00396     //got all the necessary sr info, now fill the ntuple
00397     if(!fGoodTrack) continue;
00398     
00399     fNtuple->Fill();
00400   }//end loop over events for this snarl
00401   
00402   ++fCtr;
00403 
00404   delete ntpManipulator;
00405 
00406   return result;
00407 }

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

Definition at line 1025 of file CondensedNtpModuleAtm.cxx.

Referenced by ArrivalTime_Weight().

01026 {
01027   //   cout << "in ArrivalTime_Uncertainty" << endl;
01028   
01029   Double_t sigma = 0.; // sigma is uncertainty in arrival of first p.e.
01030   
01031   // a simple monte carlo was done to calculate the uncertainty in the
01032   // arrival time of the first p.e. as a function of # of p.e.
01033   //
01034   // the only assumption made was about the shape of the pulse, a 2 ns risetime
01035   // and an 8 ns decay time of the fluor
01036   //
01037   // rlee feb 2002
01038   //
01039   
01040   if(npe<10.) sigma = 0.90+8.0/npe;
01041   else sigma = 5.0*exp(-0.5*log(npe));
01042   
01043   //  cout << sigma << endl;
01044   
01045   return sigma;
01046 }

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

Definition at line 1013 of file CondensedNtpModuleAtm.cxx.

References ArrivalTime_Uncertainty().

Referenced by GetTimeWeight().

01014 {
01015   //   cout << "in ArrivalTime_Weight" << endl;
01016   
01017   Double_t weight = 1./ArrivalTime_Uncertainty(npe);
01018   
01019   //    //cout << weight << endl;
01020   
01021   return weight*weight;
01022 }

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, Msg::kDebug, Msg::kInfo, and MSG.

00112 {
00113   // save the current working directory
00114   TDirectory* savedir  = gDirectory;  
00115   
00116     
00117   //create the new TFile for holding the electronics tree
00118   // this changes the value of gDirectory
00119   fNtpFile = new TFile(fFileName.c_str(),"RECREATE");  
00120     
00121   MSG("CondensedNtpModuleAtm", Msg::kDebug) << fFileName << " " << fTreeName 
00122                                             << " " << fNtpFile->GetName() << endl;
00123   
00124   
00125   // create TTree, these will attach themselves to the current 
00126   //working directory   
00127   fNtuple = new TTree(fTreeName.c_str(), "Analysis Tree");
00128   
00129   fNtuple->Branch("header.", "ANtpHeaderInfo", &fHeaderInfo, 64000, 2);
00130   fNtuple->Branch("beam.", "ANtpBeamInfo", &fBeamInfo, 64000, 2);
00131   fNtuple->Branch("event.", "ANtpEventInfo", &fEventInfo, 64000, 2);
00132   fNtuple->Branch("track.", "ANtpTrackInfoAtm", &fTrackInfo, 64000, 2);
00133   if(fDataType==1) fNtuple->Branch("truth.", "ANtpTruthInfoAtm", &fTruthInfo, 64000, 2);
00134   
00135   MSG("CondensedNtpModuleAtm", Msg::kDebug) << "got branches" << endl;
00136 
00137   fFailTree = new TTree("failures", "Failed Events");
00138   fFailTree->Branch("header.", "ANtpHeaderInfo", &fHeaderInfo, 64000, 2);
00139   fFailTree->Branch("multiple", &fIsMultiple, "multiple/I");
00140   fFailTree->Branch("validplane", &fValidPlaneFail, "validplane/I");
00141   fFailTree->Branch("failDeMux", &fFailDeMux, "failDeMux/I");
00142   
00143   // change back to current working directory before leaving constructor
00144   savedir->cd();   
00145   
00146   //load the vertical intensity information into some arrays
00147   //get the slant depths and the number of events for each slant depth
00148   //read in rock data
00149   // http://www.hep.umn.ed/~schubert/far/s2rock/slantdepth.data
00150   // 1-cos(zen); azi(degrees); slantdepth (gmcm**2) * cos(zen); 
00151   //    density (gm/cm**2)
00152   //  N.B., units for density are incorrectly labelled in the file
00153   ifstream densityFile(fRockMapFileName.c_str());
00154   if( !densityFile.is_open() ){
00155     cout << "could not open rock map - exit" << endl;
00156     exit(0);
00157   }
00158   
00159   int ir = 0;
00160   while( ir<1404 ){
00161     densityFile >> fPathCosZenDeg[ir] >> fPathAzimuthDeg[ir] >> fColDenCosZen[ir] >> fDensity[ir]; 
00162     ++ir;
00163   }
00164   
00165   //get the beam info tree
00166   //check to see if you need to calculate the total Protons on Target
00167   MSG("BeamTestModule", Msg::kInfo) << "beam info path = " << fBeamPath.c_str()
00168       << " tree name = " << fBeamTreeName.c_str() << endl;
00169   
00170 //   if(fDataType < 1){
00171     
00172 //     TFile *beamFile = new TFile(fBeamPath.c_str());
00173 //     fBeamTree = dynamic_cast<TTree *>(beamFile->Get(fBeamTreeName.c_str()));
00174     
00175 //     //build the index for the tree to select the utc timestamps
00176 //     fBeamTree->BuildIndex("utc");
00177 //     fBeamTree->SetBranchAddress("intensity", &fCurrentPOT);
00178 //     fBeamTree->SetBranchAddress("hornCurrent", &fHornCurrent);
00179 //     fBeamTree->SetBranchAddress("hpos", &fBeamHPos);
00180 //     fBeamTree->SetBranchAddress("vpos", &fBeamVPos);
00181 //     fBeamTree->SetBranchAddress("tgtpos", &fTargetPos);
00182 //   }
00183     
00184   return;
00185 }

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

Definition at line 1399 of file CondensedNtpModuleAtm.cxx.

Referenced by FillMCInformation().

01400 {
01401   static float Rearth = 6378140.0; /* Equatorial radius of earth */
01402   static float nuHeight = 1500.0; /* Average Neutrino production height */
01403   static float altitudeSK = -210.0; /* Altitude of the detector - MINOS is 689ft below sea level*/
01404   
01405   static float x1s = (Rearth+nuHeight)*(Rearth+nuHeight);
01406   static float x2 = Rearth+altitudeSK;
01407   static float x2s = (Rearth+altitudeSK)*(Rearth+altitudeSK);
01408   
01409   return sqrt(x1s+x2s*(cosz*cosz-1.0))-x2*cosz;
01410 }

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().

00475 {
00476   int         tmpb;  // a temp bool. See comment under DefaultConfig...
00477   int         tmpi;  // a temp int.
00478   double      tmpd;  // a temp double.
00479   const char* tmps;  // a temp string
00480   
00481   if (r.Get("FileName",           tmps)) fFileName        =  tmps;
00482   if (r.Get("RockMapFileName",    tmps)) fRockMapFileName =  tmps;
00483   if (r.Get("TreeName",           tmps)) fTreeName        =  tmps;
00484   if (r.Get("BeamPath",           tmps)) fBeamPath        = tmps;
00485   if (r.Get("BeamTreeName",       tmps)) fBeamTreeName    = tmps;
00486   if (r.Get("EndPlane",           tmpi)) fEndPlane        =  tmpi;
00487   if (r.Get("Detector",           tmpi)) fDetector        =  tmpi;
00488   if (r.Get("DataType",           tmpi)) fDataType        =  tmpi;
00489   if (r.Get("NewMC",              tmpb)) fNewMC           =  tmpb;
00490   if (r.Get("MajorRelease",       tmpi)) fMajorRelease    =  tmpi;
00491   
00492   //quiet the compiler warnings
00493   tmpb = 0;
00494   tmpd = 1.;
00495 
00496   return;
00497 }

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().

00447 {
00448   
00449   int itrue = 1;  // work around for lack of bool in registry
00450   int ifalse = 0; // work around for lack of bool in registry
00451   
00452   static Registry r;
00453   
00454   r.UnLockValues();
00455   
00456   r.Set("FileName",           "analysisTree.root");
00457   r.Set("RockMapFileName",    "density.txt");
00458   r.Set("TreeName",           "analysisTree");
00459   r.Set("BeamPath",           "/afs/fnal.gov/files/data/minos/d04/brebel/beamInfo.root");
00460   r.Set("BeamTreeName",       "beam");
00461   r.Set("EndPlane",           485);
00462   r.Set("Detector",           1);
00463   r.Set("DataType",           0);
00464   r.Set("NewMC",              itrue);
00465   r.Set("MajorRelease",       14);
00466   r.LockValues();
00467 
00468   itrue = ifalse; //quiet the compiler warnings
00469 
00470   return r;
00471 }

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, Msg::kDebug, Msg::kInfo, and MSG.

00421 {
00422 
00423   MSG("CondensedNtpModuleAtm", Msg::kDebug) << "start end job method" << endl;  
00424   
00425   //save current working directory
00426   TDirectory *savedir = gDirectory;
00427   
00428   //cd to the "directory" of the file
00429   fNtpFile->cd();
00430   
00431   MSG("CondensedNtpModuleAtm", Msg::kInfo) << fNtuple->GetEntries() << endl;
00432   
00433   fNtuple->Write();
00434   fFailTree->Write();
00435 
00436   //go back to original working directory
00437   savedir->cd();  
00438   
00439   fNtpFile->Write();
00440   fNtpFile->Close();
00441   
00442   return;
00443 }

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

Definition at line 935 of file CondensedNtpModuleAtm.cxx.

References fEventInfo, ANtpInfoObjectFiller::FillEventInformation(), fInfoFiller, Msg::kDebug, and MSG.

Referenced by Ana().

00937 {
00938   MSG("CondensedNtpModuleAtm", Msg::kDebug) << "in FillEventInformation" << endl;
00939   
00940   fInfoFiller->FillEventInformation(ntpManipulator,
00941                                     ntpEvent, 
00942                                     fEventInfo);
00943   
00944   return;
00945 }

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, Msg::kDebug, ANtpTruthInfo::leptonDCosX, ANtpTruthInfo::leptonDCosY, ANtpTruthInfo::leptonDCosZ, MSG, ANtpTruthInfo::nuEnergy, NtpMCTruth::p4neu, NtpMCStdHep::vtx, and ANtpTruthInfoAtm::weight.

Referenced by Ana().

00950 {
00951   MSG("CondensedNtpModuleAtm", Msg::kDebug) << "in FillMCInformation" << endl;
00952   
00953   const double rotMatrix[] = { 0.894507011,0,0.447053919,
00954                                0,          1,          0,
00955                               -0.447053919,0,0.894507011};
00956   
00957   double dcosx_ideal, dcosy_ideal, dcosz_ideal;
00958   
00959   fInfoFiller->FillMCTruthInformation(ntpMCTruth, fTruthInfo);
00960   
00961   if(fTruthInfo->nuEnergy>0.) fTruthInfo->baseLine = 0.001*CheapBaseline(-ntpMCTruth->p4neu[1]
00962                                                                          /TMath::Abs(fTruthInfo->nuEnergy));
00963   
00964   //get the azimuth from AstroUtill
00965   dcosx_ideal = (rotMatrix[0]*fTruthInfo->leptonDCosX
00966                  + rotMatrix[1]*fTruthInfo->leptonDCosY
00967                  + rotMatrix[2]*fTruthInfo->leptonDCosZ);
00968   dcosy_ideal = (rotMatrix[3]*fTruthInfo->leptonDCosX
00969                  + rotMatrix[4]*fTruthInfo->leptonDCosY
00970                  + rotMatrix[5]*fTruthInfo->leptonDCosZ);
00971   dcosz_ideal = (rotMatrix[6]*fTruthInfo->leptonDCosX
00972                  + rotMatrix[7]*fTruthInfo->leptonDCosY
00973                  + rotMatrix[8]*fTruthInfo->leptonDCosZ);
00974   fTruthInfo->azimuth = 180.*TMath::ATan2(-dcosx_ideal,dcosz_ideal)/TMath::Pi();
00975   if(fTruthInfo->azimuth < 0.) fTruthInfo->azimuth += 360.;
00976 
00977   //use the NtpMCStdHep object to figure out the weight for this event based on the 
00978   //fact that the new MC picks events at the surface.  we plugged the value into
00979   //the vtx Y location in the stdhep common bloc, gminos multiplies it by 10^-3
00980   //because it thought it was converting mm to m
00981   if(ntpMCStdHep){
00982     fTruthInfo->weight = 1.e3*ntpMCStdHep->vtx[1];
00983   }
00984 
00985   return;
00986 }

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, Msg::kDebug, PlaneView::kU, PlaneView::kV, ANtpTrackInfo::length, NtpSRTrack::momentum, MSG, NtpSRStrip::ndigit, NtpSRFitTrack::ndof, Munits::ns, 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().

00504 {
00505   MSG("CondensedNtpModuleAtm", Msg::kDebug) << "filling track info" << endl;
00506   
00507   //fill the basics first
00508   fInfoFiller->FillTrackInformation(ntpTrack, fTrackInfo);
00509 
00510   NtpSRStrip *ntpStrip = 0;
00511   
00512   Int_t plane = 0, strip = 0, numDoubleEndedStrips = 0, prevPlane = -1;
00513   float xPos = 0., yPos = 0.;
00514   Float_t radius = 0.;
00515   
00516   Int_t index = 0;
00517 
00518   //stuff to find the slant depth
00519   double vec1[39],vec2[36];
00520   int index1, index2, index3;
00521   
00522   for(int i=0;i<39;++i){
00523     vec1[i] = i*0.02;
00524     //     cout << vec1[i] << endl;
00525     if(i<36){
00526       vec2[i] = i*10. + 5.;
00527       //       cout << vec2[i] << endl;
00528     }
00529   }
00530 
00531   // Kashahara picks events according to where they come from
00532   //   she uses altitude, not zenith angle
00533   index1 = TMath::BinarySearch(39,vec1,1.-(-fTrackInfo->dcosYVtx))+1;
00534   index2 = TMath::BinarySearch(36,vec2,azimuth)+1;
00535   index3 = (index1)*36 +(index2);
00536   fTrackInfo->slantDepth = fDensity[index3]*0.01; //density is in 10^5g/cm^2, ie km.w.e, 
00537                                                   //0.01 converts to m.w.e 10^3m/km/(10^5g/cm^2/km.w.e)
00538   
00539   fTrackInfo->azimuth = azimuth;
00540 
00541   fTrackInfo->ra = ra;
00542   fTrackInfo->dec = dec;
00543   fTrackInfo->trackLikePlanes = ntpTrack->plane.ntrklike;
00544   Int_t nVPlanes = ntpTrack->plane.nv;
00545   Int_t nUPlanes = ntpTrack->plane.nu;
00546   fTrackInfo->invBeta = ntpTrack->time.cdtds;
00547   fTrackInfo->planesIn10Meter = 0;
00548   fTrackInfo->planesIn15Meter = 0;
00549   fTrackInfo->planesIn20Meter = 0;
00550   fTrackInfo->planesIn25Meter = 0;
00551   fTrackInfo->planesIn30Meter = 0;
00552   fTrackInfo->planesIn35Meter = 0;
00553   fTrackInfo->planesIn40Meter = 0;
00554   
00555   //check if the zenith angle indicates the right direction for the near detector. 
00556   //if not, swap the ends and adjust the azimuth
00557   if(fTrackInfo->dcosYVtx > 0. && fDetector<1){
00558     
00559     fTrackInfo->vtxX = ntpTrack->end.x; 
00560     fTrackInfo->vtxY = ntpTrack->end.y; 
00561     fTrackInfo->vtxZ = ntpTrack->end.z;                 
00562     fTrackInfo->endX = ntpTrack->vtx.x; 
00563     fTrackInfo->endY = ntpTrack->vtx.y; 
00564     fTrackInfo->endZ = ntpTrack->vtx.z; 
00565     fTrackInfo->dcosXVtx = -ntpTrack->end.dcosx;
00566     fTrackInfo->dcosYVtx = -ntpTrack->end.dcosy;
00567     fTrackInfo->dcosZVtx = -ntpTrack->end.dcosz;
00568     fTrackInfo->dcosXEnd = -ntpTrack->vtx.dcosx;
00569     fTrackInfo->dcosYEnd = -ntpTrack->vtx.dcosy;
00570     fTrackInfo->dcosZEnd = -ntpTrack->vtx.dcosz;
00571     
00572     const double rotMatrix[] = { 0.894507011,0,0.447053919,
00573                                  0,          1,          0,
00574                                  -0.447053919,0,0.894507011};
00575     
00576     double dcosx_ideal, dcosy_ideal, dcosz_ideal;
00577     
00578     //get the azimuth from AstroUtill
00579     dcosx_ideal = (rotMatrix[0]*fTrackInfo->dcosXVtx
00580                    + rotMatrix[1]*fTrackInfo->dcosYVtx
00581                    + rotMatrix[2]*fTrackInfo->dcosZVtx);
00582     dcosy_ideal = (rotMatrix[3]*fTrackInfo->dcosXVtx
00583                    + rotMatrix[4]*fTrackInfo->dcosYVtx
00584                    + rotMatrix[5]*fTrackInfo->dcosZVtx);
00585     dcosz_ideal = (rotMatrix[6]*fTrackInfo->dcosXVtx
00586                    + rotMatrix[7]*fTrackInfo->dcosYVtx
00587                    + rotMatrix[8]*fTrackInfo->dcosZVtx);
00588     
00589     fTrackInfo->azimuth = 180.*TMath::ATan2(-dcosx_ideal,dcosz_ideal)/TMath::Pi();
00590     if(fTrackInfo->azimuth < 0.) fTrackInfo->azimuth += 360.;
00591   }
00592   
00593   if(ntpTrack->fit.ndof == 0 || ntpTrack->momentum.qp== 0.){
00594     fGoodTrack = false;
00595     return;
00596   }
00597   
00598   //the value of trk.fidvtx.dz is not valid for events with the full detector
00599   //so hack up my own value of fiducialDz
00600   Float_t fullDetEndZ = 0.0594*485. + 1.;
00601   
00602   fTrackInfo->fiducialVtxDz = ntpTrack->fidvtx.dz;
00603   if(fEndPlane>248) 
00604     fTrackInfo->fiducialVtxDz = TMath::Min(fullDetEndZ - fTrackInfo->vtxZ, fTrackInfo->vtxZ); 
00605   fTrackInfo->fiducialVtxDr = ntpTrack->fidvtx.dr;
00606   
00607   fTrackInfo->fiducialEndDz = ntpTrack->fidend.dz;
00608   if(fEndPlane>248) 
00609     fTrackInfo->fiducialEndDz = TMath::Min(fullDetEndZ - fTrackInfo->endZ, fTrackInfo->endZ); 
00610   fTrackInfo->fiducialEndDr = ntpTrack->fidend.dr;
00611   
00612   fTrackInfo->uvAsymmetry = 1.;
00613   if(nVPlanes+nUPlanes>0)
00614     fTrackInfo->uvAsymmetry = TMath::Abs((1.*nUPlanes-1.*nVPlanes)/(1.*nVPlanes+1.*nUPlanes));
00615   
00616   fTrackInfo->signalUseFraction = 0.;
00617   if(fEventInfo->pulseHeight>0.) 
00618     fTrackInfo->signalUseFraction = fTrackInfo->pulseHeight/(fEventInfo->pulseHeight);
00619   
00620   fTrackInfo->planeUseFraction = 0.;
00621   if(fEventInfo->planes>0) 
00622     fTrackInfo->planeUseFraction = 1.*fTrackInfo->trackLikePlanes/(1.*fEventInfo->planes);
00623   
00624   Int_t numStrips = 0;
00625   Int_t numDigits = 0;
00626   
00627   fTrackInfo->alternateChi2 = 0.;
00628   double deltaTPos = 0.;
00629   double stripPosError = 0.041/TMath::Sqrt(12.);
00630 
00631   fTrackInfo->impactParameter = 20.0;
00632   fTrackInfo->totalStrips = 0;
00633   //does the pathlength of the track make sense?
00634   if(fTrackInfo->length<50.){
00635     
00636     
00637     //loop over all strips in the track
00638     numStrips = (int)ntpTrack->nstrip;
00639     numDoubleEndedStrips = 0;
00640     
00641     MSG("CondensedNtpModuleAtm", Msg::kDebug) << "start looping over track strips" << endl;
00642     
00643     for(Int_t ns = 0; ns < numStrips; ++ns){
00644       
00645       //get the index for this strip
00646       index = ntpTrack->stp[ns];
00647       
00648       //get the NtpSRStrip object
00649       ntpStrip = dynamic_cast<NtpSRStrip *>(stripArray->At(index));
00650       plane = (int)ntpStrip->plane;
00651       strip = (int)ntpStrip->strip;
00652       numDigits = (int)ntpStrip->ndigit;
00653       
00654       xPos = ntpTrack->stpx[ns];
00655       yPos = ntpTrack->stpy[ns];
00656      
00657       if((int)ntpStrip->planeview == PlaneView::kU)
00658         deltaTPos = ntpTrack->stpu[ns]-ntpStrip->tpos;
00659       else if((int)ntpStrip->planeview == PlaneView::kV)
00660         deltaTPos = ntpTrack->stpv[ns]-ntpStrip->tpos;
00661 
00662       MSG("CondensedNtpModuleAtm", Msg::kDebug) << ntpTrack->stpu[ns] << " " << ntpStrip->tpos
00663                                                   << " " << ntpTrack->stpv[ns] << " " << deltaTPos << endl;
00664 
00665       fTrackInfo->alternateChi2 += deltaTPos*deltaTPos/(stripPosError*stripPosError);
00666 
00667       if(numDigits >= 2){
00668         
00669         ++numDoubleEndedStrips;
00670         
00671         radius = TMath::Sqrt(yPos*yPos + xPos*xPos);
00672         if(radius<fTrackInfo->impactParameter)
00673           fTrackInfo->impactParameter = TMath::Sqrt(yPos*yPos + xPos*xPos);
00674         
00675         //get values to determine the number of plances crossed by the track 
00676         //in concentric circles in the detector from 0.3 m out to 
00677         //4.0 m.  the points come in order, so i dont have to test each one 
00678         //to make sure it is before or after the previous points
00679         if(radius>0.3 && radius<1.0 && plane != prevPlane) ++fTrackInfo->planesIn10Meter;
00680         if(radius>0.3 && radius<1.5 && plane != prevPlane) ++fTrackInfo->planesIn15Meter;
00681         if(radius>0.3 && radius<2.0 && plane != prevPlane) ++fTrackInfo->planesIn20Meter;
00682         if(radius>0.3 && radius<2.5 && plane != prevPlane) ++fTrackInfo->planesIn25Meter;
00683         if(radius>0.3 && radius<3.0 && plane != prevPlane) ++fTrackInfo->planesIn30Meter;
00684         if(radius>0.3 && radius<3.5 && plane != prevPlane) ++fTrackInfo->planesIn35Meter;
00685         if(radius>0.3 && radius<4.0 && plane != prevPlane) ++fTrackInfo->planesIn40Meter;
00686         
00687         ++fTrackInfo->totalStrips;
00688         
00689       } // end if the strip had double sided readout
00690       
00691       prevPlane = plane;
00692       
00693     } // end loop over strips in the track
00694     
00695     fTrackInfo->twoEndStripFraction = 0.;
00696     if(numStrips>0){
00697       fTrackInfo->twoEndStripFraction = (1.*numDoubleEndedStrips)/(1.*numStrips);
00698       fTrackInfo->alternateChi2 /= (1.*numStrips);
00699     }
00700   }//end if the track length is reasonable
00701   else fGoodTrack = false;
00702   
00703   MSG("CondensedNtpModuleAtm", Msg::kDebug) << "end fill track info - good track = " 
00704                                             << (int)fGoodTrack << endl;
00705   
00706   return;
00707 }

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

Definition at line 852 of file CondensedNtpModuleAtm.cxx.

References NtpSRTrack::end, fGoodTrack, fTrackInfo, Msg::kDebug, ANtpTrackInfoAtm::meanDistFromLinearFit, MSG, ANtpTrackInfoAtm::netDistFromLinearFit, Munits::ns, 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().

00854 {
00855     
00856   MSG("CondensedNtpModuleAtm", Msg::kDebug) << "in FillTrackMagneticBendingVariables" 
00857                                             << endl;
00858   
00859   NtpSRStrip *ntpStrip = 0;
00860   
00861   
00862   //find the direction cosines for the straight line between the vertex and end points
00863   //of the track
00864   Float_t deltaX = ntpTrack->end.x - ntpTrack->vtx.x;
00865   Float_t deltaY = ntpTrack->end.y - ntpTrack->vtx.y;
00866   Float_t deltaZ = ntpTrack->end.z - ntpTrack->vtx.z;
00867   Float_t deltaS = TMath::Sqrt(deltaX*deltaX + deltaY*deltaY + deltaZ*deltaZ);
00868   Float_t dcosX = deltaX/deltaS;
00869   Float_t dcosY = deltaY/deltaS;
00870   Float_t dcosZ = deltaZ/deltaS;
00871   Float_t xPos = 0.;
00872   Float_t yPos = 0.;
00873   Float_t zPos = 0.;
00874   Float_t linXPos = 0.;
00875   Float_t linYPos = 0.;
00876   
00877   //loop over the strips in the track to find the sagitta
00878   fTrackInfo->sagitta = 0.;
00879   fTrackInfo->netDistFromLinearFit = 0.;
00880   fTrackInfo->meanDistFromLinearFit = 0.;
00881   fTrackInfo->rmsDistFromLinearFit = 0.;
00882   fTrackInfo->zeroCurvatureChi2 = 0.;
00883 
00884   if(deltaZ == 0. || deltaX == 0. || deltaY == 0.){
00885     fGoodTrack = false;
00886     return;
00887   }
00888   
00889   TH1F dist("dist", "", 500, 0., 5.);
00890   
00891   Float_t curDistFromLinearFit = 0.;
00892   
00893   //loop over all strips in the track
00894   Int_t numStrips = (int)ntpTrack->nstrip;
00895   Int_t index = 0;
00896   double stripPosError = 0.041/TMath::Sqrt(12.);
00897 
00898   for(Int_t ns = 0; ns < numStrips; ++ns){
00899     
00900     //get the index for this strip
00901     index = ntpTrack->stp[ns];
00902     
00903     //get the NtpSRStrip object
00904     ntpStrip = dynamic_cast<NtpSRStrip *>(stripArray->At(index));
00905     
00906     xPos = ntpTrack->stpx[ns];
00907     yPos = ntpTrack->stpy[ns];
00908     zPos = ntpTrack->stpz[ns];
00909     
00910     linXPos = ntpTrack->vtx.x + (dcosX/dcosZ) * (zPos - ntpTrack->vtx.z);
00911     linYPos = ntpTrack->vtx.y + (dcosY/dcosZ) * (zPos - ntpTrack->vtx.z);
00912     deltaX = xPos - linXPos;
00913     deltaY = yPos - linYPos;
00914 
00915     fTrackInfo->zeroCurvatureChi2 += (deltaX*deltaX + deltaY*deltaY)/(stripPosError*stripPosError);
00916     
00917     curDistFromLinearFit = TMath::Sqrt(deltaX*deltaX + deltaY*deltaY);
00918     dist.Fill(curDistFromLinearFit);
00919     
00920     if(curDistFromLinearFit > fTrackInfo->sagitta) fTrackInfo->sagitta = curDistFromLinearFit;
00921     
00922   }
00923   
00924   if(numStrips>1) fTrackInfo->zeroCurvatureChi2 /= 1.*(numStrips-1);
00925   fTrackInfo->meanDistFromLinearFit = dist.GetMean();
00926   fTrackInfo->rmsDistFromLinearFit = dist.GetRMS();
00927   
00928   //now to figure out the net pt kick to the muon
00929   
00930   return;
00931 }

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, Msg::kDebug, SimFlag::kMC, ANtpTrackInfo::length, MSG, NtpSRStrip::ndigit, Munits::ns, 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().

00712 {    
00713   MSG("CondensedNtpModuleAtm", Msg::kDebug) << "in FillTrackTimingVariables" << endl;
00714   
00715   NtpSRStrip *ntpStrip = 0;
00716   
00717   Double_t time[2] = {0.};
00718   Double_t xPos = 0., yPos = 0., zPos = 0.;
00719   Double_t yPosition[5000] = {0.};
00720   Double_t yPosition1[5000] = {0.};
00721   Double_t weight[5000] = {0.};
00722   Double_t weight2[5000] = {0.};
00723   Double_t pathLength[5000] = {0.};
00724   Double_t weight1[5000] = {0.};
00725   Double_t timeAv[5000] = {0.};
00726   Double_t timeAv1[5000]  = {0.};
00727   Int_t arrayCtr = 0, betaCtr = 0, actr = 0, plane = 0, strip = 0, index = 0;
00728   Double_t minTime = 1.e20;     
00729   
00730   //loop over all strips in the track
00731   Int_t numStrips = (int)ntpTrack->nstrip;
00732   Int_t numDigits = 0;
00733   Double_t smear = 0.;
00734 
00735   MSG("CondensedNtpModuleAtm", Msg::kDebug) << "fill track timing" << endl;
00736   
00737   for(Int_t ns = 0; ns < numStrips; ++ns){
00738 
00739     if(fHeaderInfo->dataType == (int)SimFlag::kMC) 
00740       smear = gRandom->Gaus(0., 0.75);
00741     else smear = 0.;
00742 
00743     //get the index for this strip
00744     index = ntpTrack->stp[ns];
00745     
00746     //get the NtpSRStrip object
00747     ntpStrip = dynamic_cast<NtpSRStrip *>(stripArray->At(index));
00748     plane = (int)ntpStrip->plane;
00749     strip = (int)ntpStrip->strip;
00750     numDigits = (int)ntpStrip->ndigit;
00751     
00752     xPos = ntpTrack->stpx[ns];
00753     yPos = ntpTrack->stpy[ns];
00754     zPos = ntpTrack->stpz[ns];
00755     
00756     if(numDigits == 2){
00757       
00758       time[0] = ntpTrack->stpt0[ns];
00759       time[1] = ntpTrack->stpt1[ns];
00760       
00761       if(time[0] < minTime) minTime = time[0];
00762       if(time[1] < minTime) minTime = time[1];
00763       
00764       //have a point for each end of the strip
00765       timeAv[arrayCtr] = 1.e9*time[0] + smear;
00766       timeAv[arrayCtr+1] = 1.e9*time[1] + smear;
00767       pathLength[arrayCtr] = fTrackInfo->length - ntpTrack->stpds[ns];
00768       pathLength[arrayCtr+1] = fTrackInfo->length - ntpTrack->stpds[ns];
00769       
00770       //stuff to get timing information
00771       yPosition[arrayCtr] = yPos;
00772       weight[arrayCtr] = GetTimeWeight(ntpStrip->ph0.sigcor);
00773       if (weight[arrayCtr]>0.4) weight[arrayCtr]=0.4;
00774       yPosition[arrayCtr+1] = yPos;
00775       weight[arrayCtr+1] = GetTimeWeight(ntpStrip->ph1.sigcor);
00776       if (weight[arrayCtr+1]>0.4) weight[arrayCtr+1]=0.4;
00777       MSG("CondensedNtpModuleAtm", Msg::kDebug) << weight[arrayCtr] << " " 
00778                                                 << weight[arrayCtr+1] << endl;
00779       arrayCtr += 2;
00780       
00781       
00782     } // end if the strip had double sided readout
00783   } // end loop over strips in the track
00784 
00785   //subtract the minimum time and copy the weights into another array to do the fit to 
00786   //a flat line later
00787   
00788   for(Int_t ti = 0; ti< arrayCtr; ++ti){
00789     timeAv[ti] -= 1.e9*minTime;
00790     weight2[ti] = weight[ti];
00791   }
00792   
00793   //cout << "got strip info" << endl;
00794   
00795   //array for fit parameters
00796   Double_t param[] = {-10000., -10000.};
00797   Double_t eparam[] = {-10000., -10000.};
00798   
00799   //find the inverse beta value
00800   betaCtr = arrayCtr;
00801   FitMinimizingResidual(betaCtr, pathLength, timeAv, 
00802                         weight, param, eparam);
00803   fTrackInfo->invBeta = param[1]*0.3;
00804   if(betaCtr > 2) fTrackInfo->invBetaChi2 = eparam[0] /(1.*(betaCtr-2));
00805   
00806   //find the chi^2 for a fit to a straight line so you 
00807   //can see if the sloped fit is really better
00808   FitMinimizingResidual(betaCtr, pathLength, timeAv, 
00809                         weight2, param, eparam, false);
00810   if(betaCtr > 1) fTrackInfo->flatLineChi2 = eparam[0] / (1.*betaCtr-1);
00811 
00812   //now loop over the points in the array and fill new arrays 
00813   //containing only those points that have a nonzero weight
00814   actr = 0;
00815   MSG("CondensedNtpModuleAtm", Msg::kDebug) << fHeaderInfo->run << " " 
00816                                             << fHeaderInfo->subRun
00817                                             << " " << fHeaderInfo->snarl 
00818                                             << " " << arrayCtr << endl;
00819   for(Int_t ii = 0; ii < arrayCtr; ++ii){
00820     MSG("CondensedNtpModuleAtm", Msg::kDebug) << " " << ii << "/" << arrayCtr << " " 
00821                                               << yPosition[ii] << " " << weight[ii] 
00822                                               << " " << timeAv[ii] << " " 
00823                                               << actr << endl;
00824     if(weight[ii]>0.){
00825       yPosition1[actr] = yPosition[ii];
00826       weight1[actr] = weight[ii];
00827       timeAv1[actr] = timeAv[ii];
00828       ++actr;
00829     }
00830   }//end loop to select good points along the track
00831   
00832   //get the fit for the time in the y direction
00833   FitMinimizingResidual(actr, yPosition1, timeAv1, weight1, param, eparam);
00834   fTrackInfo->timeSlopeY = param[1];
00835   fTrackInfo->timeSlopeYChi2 = eparam[0];
00836   
00837   if(actr > 2) fTrackInfo->timeSlopeYChi2 /= 1.*(actr-2);  
00838   
00839   MSG("CondensedNtpModuleAtm", Msg::kDebug) << "actr = " << actr 
00840                                             << "/" << arrayCtr 
00841                                             << " goodtrack = " 
00842                                             << (int)fGoodTrack << endl;
00843   
00844   if(1.*actr < 0.5*arrayCtr) fGoodTrack = false;
00845   
00846   return;
00847 }

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.

01281 {
01282   
01283   double vtxX = x[0];
01284   double vtxY = y[0];
01285   double endX = x[nPoints-1];
01286   double endY = y[nPoints-1];
01287   
01288   a2 = endY-vtxY;
01289   if(TMath::Abs(vtxX-endX)>0.) a2 /= (endX-vtxX);
01290   else{
01291     cout << "cannot find slope for n = " << nPoints << " points " << endl;
01292     for(Int_t i = 0; i<nPoints; ++i) 
01293       cout << x[i] << " " << y[i] << " " << weight[i] << endl;
01294   }
01295   
01296   a1 = endY-a2*endX;
01297   
01298   //find the chi^2 value for the fit
01299   chiSq = 0.;
01300   Double_t chi = 0.;
01301   for(Int_t j = 0; j<nPoints; j++){
01302     
01303     chi = (y[j] - (a1 + a2*x[j]))/weight[j];
01304     
01305     //     cout << y[j] << " " << a1+a2*x[j] << " " << chiSq << endl;
01306     
01307     chiSq += chi*chi;
01308     
01309   }  
01310   
01311   //   cout << "chi^2 = " << chiSq << " " << a1 << " " << a2 << endl;
01312   
01313   //get the chiSq per degree of freedom
01314   
01315   if(nPoints>2) chiSq /= 1.*(nPoints-2);
01316   
01317   return;
01318 }

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.

01183 {
01184   
01185   //use method of determinants to do the least squares fit
01186   //this is from bevington chapter 6
01187   //
01188   //f_1(x) = 1
01189   //f_2(x) = x
01190   
01191   //need to find the sums
01192   //SIGMA y_i*f_1(x_i)/weight_i
01193   //SIGMA y_i*f_2(x_i)/weight_i
01194   //SIGMA f_1(x_i)*f_1(x_i)/weight_i
01195   //SIGMA f_1(x_i)*f_2(x_i)/weight_i
01196   //SIGMA f_2(x_i)*f_1(x_i)/weight_i
01197   //SIGMA f_2(x_i)*f_2(x_i)/weight_i
01198   //
01199   //some of these are the same, ie f_2*f_1 = f_1*f_2 so only need
01200   //3 instead of 4 variables
01201   
01202   Double_t yf1divWeight = 0.;
01203   Double_t yf2divWeight = 0.;
01204   Double_t f1f1divWeight = 0.;
01205   Double_t f1f2divWeight = 0.;
01206   Double_t f2f2divWeight = 0.;
01207   
01208   for(Int_t i=0; i<nPoints; i++){
01209     //     cout << y[i] << " " << x[i] << " " << weight[i] << endl;
01210     
01211     yf1divWeight += y[i]/(weight[i]*weight[i]);
01212     yf2divWeight += (y[i]*x[i])/(weight[i]*weight[i]);
01213     f1f1divWeight += 1./(weight[i]*weight[i]);
01214     f1f2divWeight += x[i]/(weight[i]*weight[i]);
01215     f2f2divWeight += (x[i]*x[i])/(weight[i]*weight[i]);
01216   }
01217   
01218   //make a bunch of matricies to hold these values
01219   TMatrix deltaMatrix(2,2);
01220   TMatrix a1Matrix(2,2);
01221   TMatrix a2Matrix(2,2);
01222   
01223   //fill the matricies
01224   deltaMatrix(0,0) = f1f1divWeight;
01225   deltaMatrix(0,1) = f1f2divWeight;
01226   deltaMatrix(1,0) = f1f2divWeight;
01227   deltaMatrix(1,1) = f2f2divWeight;
01228   
01229   a1Matrix(0,0) = yf1divWeight;
01230   a1Matrix(0,1) = f1f2divWeight;
01231   a1Matrix(1,0) = yf2divWeight;
01232   a1Matrix(1,1) = f2f2divWeight;
01233   
01234   a2Matrix(0,0) = f1f1divWeight;
01235   a2Matrix(0,1) = yf1divWeight;
01236   a2Matrix(1,0) = f1f2divWeight;
01237   a2Matrix(1,1) = yf2divWeight;
01238   
01239   a1 = -10000.;
01240   a2 = -10000.;
01241   Double_t delta = deltaMatrix.Determinant();
01242   
01243   //   cout << "delta = " << delta << endl;
01244   
01245   if(delta != 0.){
01246     a1 = a1Matrix.Determinant()/delta;
01247     a2 = a2Matrix.Determinant()/delta;
01248   }
01249   
01250   //find the chi^2 value for the fit
01251   chiSq = 0.;
01252   Double_t chi = 0.;
01253   for(Int_t j = 0; j<nPoints; j++){
01254     
01255     chi = (y[j] - (a1 + a2*x[j]))/weight[j];
01256     
01257     //     cout << chi*weight[j] << endl;
01258     
01259     chiSq += chi*chi;
01260     
01261   }  
01262   
01263   //   cout << "chi^2 = " << chiSq << " " << a1 << " " << a2 << endl;
01264   
01265   //get the chiSq per degree of freedom
01266   
01267   chiSq /= (1.*nPoints - 2.);
01268   
01269   return;
01270 }

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

Definition at line 1366 of file CondensedNtpModuleAtm.cxx.

01369 {
01370   
01371   //   cout << "rms" << endl;
01372   
01373   double sumXXW = 0.;
01374   double sumXW = 0.;
01375   double sumW = 0.;
01376   
01377   for(Int_t i = 0; i<nPoints; ++i){
01378     sumXXW += x[i]*x[i]*weight[i];
01379     sumXW += x[i]*weight[i];
01380     sumW += weight[i];
01381     //     cout << sumXW << " " << sumXXW << " " << sumW << " " 
01382     //   << x[i] << " " << weight[i] << endl;
01383   }
01384   
01385   if(sumW > 0. && nPoints>1){
01386     mean = sumXW/sumW;
01387     if(sumXXW >=0.) rms =  TMath::Sqrt((sumXXW/sumW - mean*mean)
01388                                        *(nPoints*1./(1.*nPoints-1.)));
01389   }
01390   else{
01391     rms = 0.;
01392     mean = x[0];
01393   }
01394   return;
01395 }

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.

01331 {
01332   
01333   double vtxX = x[0];
01334   double vtxY = y[0];
01335   double endX = x[nPoints-1];
01336   double endY = y[nPoints-1];
01337   
01338   a2 = endY-vtxY;
01339   if(TMath::Abs(vtxX-endX)>0.) a2 /= (endX-vtxX);
01340   else{
01341     cout << "cannot find slope for n = " << nPoints << " points " << endl;
01342     for(Int_t i = 0; i<nPoints; ++i) cout << x[i] << " " << y[i] << endl;
01343   }
01344   
01345   a1 = endY-a2*endX;
01346   
01347   //find the chi^2 value for the fit
01348   deviation = 0.;
01349   Double_t dev = 0.;
01350   for(Int_t j = 0; j<nPoints; j++){
01351     
01352     dev = (y[j] - (a1 + a2*x[j]));
01353     
01354     //     cout << y[j] << " " << a1+a2*x[j] << " " << chiSq << endl;
01355     
01356     deviation += dev*dev;
01357     
01358   }  
01359   
01360   return;
01361 }

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().

01123 {
01124   //   cout << "in LinearFit_MinimizeResid" << endl;
01125   
01126   Int_t nremove = 0;
01127   Double_t maxresid = 11.;
01128   Double_t resid = 0.;
01129   Double_t parm[2] = {0.};
01130   Double_t eparm[2] = {0.};
01131   Double_t timefitchi2 = 0.;
01132   while(maxresid>10.
01133         &&n-nremove>4
01134         &&(1.*nremove)/(1.*n)<0.2){
01135     
01136     timefitchi2=0.;
01137     if(findSlope) LinearFit_Weighted(n,x,y,w,parm,eparm);
01138     else WeightedAverage(n,y,w,parm,eparm);
01139     maxresid = 0.;
01140     Int_t imaxresid=0;
01141     for (Int_t i=0; i<n; i++) {
01142       if (w[i]>0.) {
01143         if(findSlope) resid = (parm[1]*x[i]+parm[0]-y[i]);
01144         else resid = parm[0] - y[i];
01145         timefitchi2 += w[i]*resid*resid;
01146         if (TMath::Abs(resid)>maxresid) {
01147           maxresid = TMath::Abs(resid);
01148           imaxresid = i;
01149         }
01150       }
01151     }
01152     if (maxresid>10.) {
01153       w[imaxresid] = 0.;
01154       nremove++;
01155     }
01156   }//end loop to minimize residuals
01157   
01158   //get the chi^2
01159   timefitchi2 = 0.;
01160   for (Int_t i=0; i<n; i++) {
01161     if (w[i]>0.) {
01162       if(findSlope) resid = (parm[1]*x[i]+parm[0]-y[i]);
01163       else resid = parm[0] - y[i];
01164       timefitchi2 += w[i]*resid*resid;
01165     }
01166   }
01167   
01168   
01169   n -= nremove;
01170   param[0] = parm[0];
01171   param[1] = parm[1];
01172   eparam[0] = timefitchi2;
01173   eparam[1] = eparm[1];
01174   
01175   return;
01176 }

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

Definition at line 1000 of file CondensedNtpModuleAtm.cxx.

References ArrivalTime_Weight().

Referenced by FillTrackTimingVariables().

01001 {
01002   //   cout << "in GetTimeWeight" << endl;
01003   
01004   Double_t npe = ph/60.;
01005   if (npe<1.) npe=1.;
01006   
01007   // cout << ArrivalTime_Weight(npe) << endl;
01008   
01009   return ArrivalTime_Weight(npe);
01010 }

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::kInfo, and MSG.

00411 {
00412   MSG("JobC", Msg::kInfo) 
00413     << "NearElectronicsCheck::Help\n"
00414     <<"CondensedNtpModuleAtm is a module which demultiplexes events "
00415     <<"in the far detector."
00416     << endl;
00417 }

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().

01055 {
01056   //   cout << "in LinearFit_Weighted with errors" << endl;
01057   
01058   Double_t sumx=0.;
01059   Double_t sumx2=0.;
01060   Double_t sumy=0.;
01061   Double_t sumy2=0.;
01062   Double_t sumxy=0.;
01063   Double_t sumw=0.;
01064   
01065   for(Int_t i=0; i<n; ++i){
01066     sumx += x[i]*w[i];
01067     sumx2 += x[i]*x[i]*w[i];
01068     sumy += y[i]*w[i]; 
01069     sumy2 += y[i]*y[i]*w[i];
01070     sumxy += x[i]*y[i]*w[i];
01071     sumw += w[i];
01072     
01073     //     cout << x[i] << " " << y[i] << " " << w[i] << endl;
01074   }
01075   
01076   //   cout << sumx2 << " " << sumw << " " << sumx << " " << sumxy << " " << sumy
01077   //        << endl;
01078   
01079   if(sumx2*sumw-sumx*sumx>0.&&sumx2<1.e5){
01080     parm[0] = (sumy*sumx2-sumx*sumxy)/(sumx2*sumw-sumx*sumx);
01081     parm[1] = (sumxy*sumw-sumx*sumy)/(sumx2*sumw-sumx*sumx);
01082     
01083     eparm[0] = TMath::Sqrt(sumx2*(sumx2*sumw-sumx*sumx))/(sumx2*sumw-sumx*sumx);
01084     eparm[1] = TMath::Sqrt(sumx2*sumw-sumx*sumx)/(sumx2*sumw-sumx*sumx);
01085   }
01086   
01087   return;
01088 }

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().

00990 {
00991   fEventInfo->Reset();
00992   fTrackInfo->Reset();
00993   fTruthInfo->Reset();
00994   fGoodTrack = false;
00995   
00996   return;
00997 }

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().

01095 {
01096   //   cout << "in LinearFit_Weighted with errors" << endl;
01097 
01098   double sumWeights = 0.;
01099   double sumXWeights = 0.;
01100   //find the weighted average of the points
01101   for(int i = 0; i < n; ++i){
01102     sumXWeights += x[i]*w[i];
01103     sumWeights += w[i];
01104   }
01105   
01106   if(sumWeights > 0.) parm[0] = sumXWeights/sumWeights;
01107   else parm[0] = -9999.;
01108   parm[1] = 0.;
01109   
01110   eparm[0] = 0.;
01111   eparm[1] = 0.;
01112 
01113   
01114   return;
01115 }


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(), BeginJob(), and ~CondensedNtpModuleAtm().

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.

Referenced by Ana().

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 Ana(), 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(), ResetTreeVariables(), and ~CondensedNtpModuleAtm().

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(), FillTrackTimingVariables(), and ~CondensedNtpModuleAtm().

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 Ana(), and 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(), ResetTreeVariables(), and ~CondensedNtpModuleAtm().

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(), ResetTreeVariables(), and ~CondensedNtpModuleAtm().

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:
Generated on Thu Jul 10 22:52:36 2014 for loon by  doxygen 1.4.7