MeuAnalysis Class Reference

#include <MeuAnalysis.h>

List of all members.

Public Member Functions
	MeuAnalysis ()
	~MeuAnalysis ()
void	Test ()
void	TestEventLoop ()
void	WriteOutHistos () const
void	BasicPlots ()
void	BasicReco ()
void	MakeSummaryTree ()
void	N_1Plots ()
void	SnarlList ()
void	SpillPlots ()
void	MakeSummaryTreeWithNtpStOneSnarl (const NtpStRecord *pntp, const NtpBDLiteRecord *pntpBD) const
void	StoreOrFinishSummaryTree (MeuSummaryWriter *psOut, MeuCounter *pcnt, Bool_t finish) const
Bool_t	ObjectExistsInFile (TFile *file, std::string objectName) const
std::string	GetFirstFileName (std::string wildcardString) const
Static Public Member Functions
void	InputFileName (std::string f)
void	LoadTrees (Bool_t loadTrees)
void	UseAtNu ()
void	RecalibrateSigLin (Bool_t recalibrate=true)
void	RecalibrateSigCor (Bool_t recalibrate=true)
void	RecalibratePE (Bool_t recalibrate=true)
Private Member Functions
Float_t	GetTemperature (const MeuSummary &s) const
void	InitialiseLoopVariables ()
void	MakeChain ()
void	MakeSummaryTreeWithAtNu ()
void	MakeSummaryTreeWithNtpSt ()
std::vector< std::string >	MakeFileList ()
void	SetRootFileObjects ()
TFile *	OpenFile (Int_t runNumber, std::string prefix) const
ofstream *	OpenTxtFile (Int_t runNumber, std::string prefix) const
void	SetLoopVariables (Int_t event, Int_t printMode=1)
Bool_t	PlaneIsWithinTrack (const MeuSummary &s, Int_t plane, Int_t vtxPl, Int_t endPl) const
void	Recalibrate (const NtpStRecord &ntp, MeuSummary &s, std::map< Int_t, MeuHitInfo > &plInfo, const NtpSREvent &evt) const
void	RecalibrateAtNu (AtmosEvent &ev, MeuSummary &s, std::map< Int_t, MeuHitInfo > &plInfo) const
void	CheckTrackGaps2 (const NtpStRecord &ntp) const
Bool_t	FilterLI (const NtpStRecord &ntp, const MeuSummary &s) const
Bool_t	FilterLI (const AtmosEvent &ev) const
Private Attributes
TChain *	fChain
TChain *	fChainBD
AtmosEvent *	fev
TFile *	fOutFile
NtpStRecord *	fRec
NtpBDLiteRecord *	fRecBD
Bool_t	fUseAtNu
Float_t	fEntries
Float_t	fEntriesBD
Int_t	fRun
Int_t	fSubrun
Static Private Attributes
std::string	fInputFileName = ""
Bool_t	fLoadTrees = false
Bool_t	fRecalibrateSigLin = false
Bool_t	fRecalibrateSigCor = false
Bool_t	fRecalibratePE = false

---

Constructor & Destructor Documentation

MeuAnalysis::MeuAnalysis (   )

Definition at line 71 of file MeuAnalysis.cxx. References fChain, RecCandHeader::GetEvent(), RecRecordImp< T >::GetHeader(), RecDataHeader::GetRun(), RecDataHeader::GetSubRun(), and MSG. { MSG("MeuAnalysis",Msg::kInfo) <<"Running MeuAnalysis Constructor..."<<endl; fChain=0; fChainBD=0; fev=0;//atnu one fOutFile=0; fRec=0;//sr fRecBD=0;//sr beam data fUseAtNu=false; fEntries=0; fEntriesBD=0; fRun=100;//default fSubrun=0; if (fLoadTrees) { MSG("MeuAnalysis",Msg::kInfo)<<"Loading trees"<<endl; this->MakeChain(); this->SetRootFileObjects(); //this->InitialiseHitInfoVariables(); //get basic run number stuff if (fRec && fChain){//sr //get the first event fChain->GetEvent(0); NtpStRecord& ntp=(*fRec); const RecCandHeader& rec=ntp.GetHeader(); MSG("MeuAnalysis",Msg::kInfo) <<"Found: run="<<rec.GetRun() <<", subrun="<<rec.GetSubRun()<<endl; fRun=rec.GetRun(); fSubrun=rec.GetSubRun(); } else if (fev && fChain){//atnu //get the first event fChain->GetEvent(0); MSG("MeuAnalysis",Msg::kInfo) <<"Found: run="<<fev->Run <<", subrun="<<fev->SubRun<<endl; fRun=fev->Run; fSubrun=fev->SubRun; } else cout<<"Ahhh, no root files"<<endl; } else { MSG("MeuAnalysis",Msg::kInfo) <<"Not loading trees in constructor"<<endl; } //setup some root stuff //include the under and overflow counts gStyle->SetOptStat(1111111); gStyle->SetOptFit(1111); gStyle->SetPalette(1); MSG("MeuAnalysis",Msg::kInfo) <<"Finished MeuAnalysis Constructor"<<endl; }

MeuAnalysis::~MeuAnalysis (   )

Definition at line 136 of file MeuAnalysis.cxx. References fOutFile, and MSG. { MSG("MeuAnalysis",Msg::kInfo) <<"Running MeuAnalysis Destructor..."<<endl; if (fOutFile){ //this makes histos disappear from the canvases //so do it at the very end fOutFile->Close(); } MSG("MeuAnalysis",Msg::kDebug) <<"Finished MeuAnalysis Destructor"<<endl; }

---

Member Function Documentation

void MeuAnalysis::BasicPlots (   )

Definition at line 1607 of file MeuAnalysis.cxx. References NtpSRFiducial::dr, NtpSRTrack::ds, NtpSRFiducial::dz, NtpSRTrack::end, NtpStRecord::evthdr, NtpSRTrack::fidend, FilterLI(), Form(), fOutFile, fRec, InitialiseLoopVariables(), MAXMSG, MSG, NtpSREventSummary::ndigit, NtpSRTrack::nstrip, NtpSREventSummary::ntrack, OpenFile(), NtpSRStrip::ph1, MeuHitInfo::Plane, NtpSRStrip::plane, NtpSRVertex::plane, NtpSRStrip::planeview, NtpSRTrack::range, MeuSummary::RFid, s(), SetLoopVariables(), NtpSRPulseHeight::sigcor, MeuHitInfo::SigCor2, NtpSRTrack::stp, NtpStRecord::stp, NtpSRTrack::stpx, NtpSRTrack::stpy, NtpSRTrack::stpz, MeuHitInfo::Strip, NtpSRStrip::strip, MeuHitInfo::TPos, NtpSRStrip::tpos, NtpStRecord::trk, MeuHitInfo::View, NtpSRTrack::vtx, MeuHitInfo::X, NtpSRVertex::x, MeuHitInfo::Y, NtpSRVertex::y, MeuHitInfo::Z, NtpSRStrip::z, and NtpSRVertex::z. { MSG("MeuAnalysis",Msg::kInfo) <<" ** Running BasicPlots method... **"<<endl; //open the output file for the histograms fOutFile=this->OpenFile(100,"BasicPlots"); TH1F* hNTrack=new TH1F("hNTrack","hNTrack",10,-1,9); hNTrack->SetTitle("NTrack"); hNTrack->GetXaxis()->SetTitle("NTrack"); hNTrack->GetXaxis()->CenterTitle(); hNTrack->GetYaxis()->SetTitle(""); hNTrack->GetYaxis()->CenterTitle(); hNTrack->SetFillColor(0); //hNTrack->SetBit(TH1::kCanRebin); TH1F* hStPerTrk=new TH1F("hStPerTrk","hStPerTrk",201,-1,200); hStPerTrk->SetTitle("StPerTrk"); hStPerTrk->GetXaxis()->SetTitle("StPerTrk"); hStPerTrk->GetXaxis()->CenterTitle(); hStPerTrk->GetYaxis()->SetTitle(""); hStPerTrk->GetYaxis()->CenterTitle(); hStPerTrk->SetFillColor(0); //hStPerTrk->SetBit(TH1::kCanRebin); TH1F* hTPos=new TH1F("hTPos","hTPos",200,-4,4); hTPos->SetTitle("TPos"); hTPos->GetXaxis()->SetTitle("TPos"); hTPos->GetXaxis()->CenterTitle(); hTPos->GetYaxis()->SetTitle(""); hTPos->GetYaxis()->CenterTitle(); hTPos->SetFillColor(0); //hTPos->SetBit(TH1::kCanRebin); TH1F* hSigCor=new TH1F("hSigCor","hSigCor",3000,-1,15000); hSigCor->SetTitle("SigCor"); hSigCor->GetXaxis()->SetTitle("SigCor"); hSigCor->GetXaxis()->CenterTitle(); hSigCor->GetYaxis()->SetTitle(""); hSigCor->GetYaxis()->CenterTitle(); hSigCor->SetFillColor(0); //hSigCor->SetBit(TH1::kCanRebin); TH1F* hView=new TH1F("hView","hView",6,-1,5); hView->SetTitle("View"); hView->GetXaxis()->SetTitle("View"); hView->GetXaxis()->CenterTitle(); hView->GetYaxis()->SetTitle(""); hView->GetYaxis()->CenterTitle(); hView->SetFillColor(0); //hView->SetBit(TH1::kCanRebin); TH1F* hPlane=new TH1F("hPlane","hPlane",301,-1,300); hPlane->SetTitle("Plane"); hPlane->GetXaxis()->SetTitle("Plane"); hPlane->GetXaxis()->CenterTitle(); hPlane->GetYaxis()->SetTitle(""); hPlane->GetYaxis()->CenterTitle(); hPlane->SetFillColor(0); //hPlane->SetBit(TH1::kCanRebin); TH1F* hStrip=new TH1F("hStrip","hStrip",101,-1,100); hStrip->SetTitle("Strip"); hStrip->GetXaxis()->SetTitle("Strip"); hStrip->GetXaxis()->CenterTitle(); hStrip->GetYaxis()->SetTitle(""); hStrip->GetYaxis()->CenterTitle(); hStrip->SetFillColor(0); //hStrip->SetBit(TH1::kCanRebin); TH1F* hZ=new TH1F("hZ","hZ",100,-1,30); hZ->SetTitle("Z"); hZ->GetXaxis()->SetTitle("Z"); hZ->GetXaxis()->CenterTitle(); hZ->GetYaxis()->SetTitle(""); hZ->GetYaxis()->CenterTitle(); hZ->SetFillColor(0); //hZ->SetBit(TH1::kCanRebin); TH1F* hR=new TH1F("hR","hR",100,-1,10); hR->SetTitle("R"); hR->GetXaxis()->SetTitle("R"); hR->GetXaxis()->CenterTitle(); hR->GetYaxis()->SetTitle(""); hR->GetYaxis()->CenterTitle(); hR->SetFillColor(0); //hR->SetBit(TH1::kCanRebin); TH1F* hMaxR=new TH1F("hMaxR","hMaxR",100,-1,10); hMaxR->SetTitle("MaxR"); hMaxR->GetXaxis()->SetTitle("MaxR"); hMaxR->GetXaxis()->CenterTitle(); hMaxR->GetYaxis()->SetTitle(""); hMaxR->GetYaxis()->CenterTitle(); hMaxR->SetFillColor(0); //hMaxR->SetBit(TH1::kCanRebin); TH1F* hMinR=new TH1F("hMinR","hMinR",100,-1,10); hMinR->SetTitle("MinR"); hMinR->GetXaxis()->SetTitle("MinR"); hMinR->GetXaxis()->CenterTitle(); hMinR->GetYaxis()->SetTitle(""); hMinR->GetYaxis()->CenterTitle(); hMinR->SetFillColor(0); //hMinR->SetBit(TH1::kCanRebin); TH2F* hTPosVsPlane=new TH2F("hTPosVsPlane","hTPosVsPlane", 301,-1,300,100,-3,3); hTPosVsPlane->SetTitle("TPos vs Plane"); hTPosVsPlane->GetXaxis()->SetTitle("Plane"); hTPosVsPlane->GetXaxis()->CenterTitle(); hTPosVsPlane->GetYaxis()->SetTitle("TPos"); hTPosVsPlane->GetYaxis()->CenterTitle(); hTPosVsPlane->SetFillColor(0); //hTPosVsPlane->SetBit(TH1::kCanRebin); TH2F* hYvsX=new TH2F("hYvsX","hYvsX", 100,-5,5,100,-5,5); hYvsX->SetTitle("Y vs X"); hYvsX->GetXaxis()->SetTitle("X"); hYvsX->GetXaxis()->CenterTitle(); hYvsX->GetYaxis()->SetTitle("Y"); hYvsX->GetYaxis()->CenterTitle(); hYvsX->SetFillColor(0); //hYvsX->SetBit(TH1::kCanRebin); TH2F* hYvsXVtx=new TH2F("hYvsXVtx","hYvsXVtx", 100,-5,5,100,-5,5); hYvsXVtx->SetTitle("Y vs X"); hYvsXVtx->GetXaxis()->SetTitle("X"); hYvsXVtx->GetXaxis()->CenterTitle(); hYvsXVtx->GetYaxis()->SetTitle("Y"); hYvsXVtx->GetYaxis()->CenterTitle(); hYvsXVtx->SetFillColor(0); //hYvsXVtx->SetBit(TH1::kCanRebin); TH2F* hYvsXEnd=new TH2F("hYvsXEnd","hYvsXEnd", 100,-5,5,100,-5,5); hYvsXEnd->SetTitle("Y vs X"); hYvsXEnd->GetXaxis()->SetTitle("X"); hYvsXEnd->GetXaxis()->CenterTitle(); hYvsXEnd->GetYaxis()->SetTitle("Y"); hYvsXEnd->GetYaxis()->CenterTitle(); hYvsXEnd->SetFillColor(0); //hYvsXEnd->SetBit(TH1::kCanRebin); TH2F* hYvsXVtxFid=new TH2F("hYvsXVtxFid","hYvsXVtxFid", 100,-5,5,100,-5,5); hYvsXVtxFid->SetTitle("Y vs X"); hYvsXVtxFid->GetXaxis()->SetTitle("X"); hYvsXVtxFid->GetXaxis()->CenterTitle(); hYvsXVtxFid->GetYaxis()->SetTitle("Y"); hYvsXVtxFid->GetYaxis()->CenterTitle(); hYvsXVtxFid->SetFillColor(0); //hYvsXVtxFid->SetBit(TH1::kCanRebin); TH2F* hYvsXEndFid=new TH2F("hYvsXEndFid","hYvsXEndFid", 100,-5,5,100,-5,5); hYvsXEndFid->SetTitle("Y vs X"); hYvsXEndFid->GetXaxis()->SetTitle("X"); hYvsXEndFid->GetXaxis()->CenterTitle(); hYvsXEndFid->GetYaxis()->SetTitle("Y"); hYvsXEndFid->GetYaxis()->CenterTitle(); hYvsXEndFid->SetFillColor(0); //hYvsXEndFid->SetBit(TH1::kCanRebin); TH2F* hYvsZ=new TH2F("hYvsZ","hYvsZ",400,-1,30,100,-5,5); hYvsZ->SetTitle("Y vs Z"); hYvsZ->GetXaxis()->SetTitle("Z"); hYvsZ->GetXaxis()->CenterTitle(); hYvsZ->GetYaxis()->SetTitle("Y"); hYvsZ->GetYaxis()->CenterTitle(); hYvsZ->SetFillColor(0); //hYvsZ->SetBit(TH1::kCanRebin); TH2F* hXvsZ=new TH2F("hXvsZ","hXvsZ",400,-1,30,100,-5,5); hXvsZ->SetTitle("X vs Z"); hXvsZ->GetXaxis()->SetTitle("Z"); hXvsZ->GetXaxis()->CenterTitle(); hXvsZ->GetYaxis()->SetTitle("Y"); hXvsZ->GetYaxis()->CenterTitle(); hXvsZ->SetFillColor(0); //hXvsZ->SetBit(TH1::kCanRebin); TH2F* hYvsPlane=new TH2F("hYvsPlane","hYvsPlane",301,-1,300,100,-5,5); hYvsPlane->SetTitle("Y vs Z"); hYvsPlane->GetXaxis()->SetTitle("Z"); hYvsPlane->GetXaxis()->CenterTitle(); hYvsPlane->GetYaxis()->SetTitle("Y"); hYvsPlane->GetYaxis()->CenterTitle(); hYvsPlane->SetFillColor(0); //hYvsPlane->SetBit(TH1::kCanRebin); TH2F* hXvsPlane=new TH2F("hXvsPlane","hXvsPlane",301,-1,300,100,-5,5); hXvsPlane->SetTitle("X vs Z"); hXvsPlane->GetXaxis()->SetTitle("Z"); hXvsPlane->GetXaxis()->CenterTitle(); hXvsPlane->GetYaxis()->SetTitle("Y"); hXvsPlane->GetYaxis()->CenterTitle(); hXvsPlane->SetFillColor(0); //hXvsPlane->SetBit(TH1::kCanRebin); vector<TH1F*> vhStrips(282); vector<TH1F*> vhTPos(282); for (Int_t i=0;i<282;i++) { string sName="hStripsPlane"; string sPl=Form("%d",i); sName+=sPl; //cout<<"name="<<sName<<endl; TH1F* hStrip=new TH1F(sName.c_str(),sName.c_str(), 101,-1,100); hStrip->SetTitle(sName.c_str()); hStrip->GetXaxis()->SetTitle("Strip"); hStrip->GetXaxis()->CenterTitle(); hStrip->GetYaxis()->SetTitle(""); hStrip->GetYaxis()->CenterTitle(); hStrip->SetFillColor(0); //hStrip->SetBit(TH1::kCanRebin); vhStrips[i]=hStrip; //strips are 0-95(96), 4-67(64), 0-63(64) sName="hTPosPlane"; sName+=sPl; TH1F* hTPos=new TH1F(sName.c_str(),sName.c_str(),200,-4,4); hTPos->SetTitle(sName.c_str()); hTPos->GetXaxis()->SetTitle("TPos"); hTPos->GetXaxis()->CenterTitle(); hTPos->GetYaxis()->SetTitle(""); hTPos->GetYaxis()->CenterTitle(); hTPos->SetFillColor(0); //hTPos->SetBit(TH1::kCanRebin); vhTPos[i]=hTPos; } vector<TH1F*> vhStripsTrk(282); vector<TH1F*> vhTPosTrk(282); for (Int_t i=0;i<282;i++) { string sName="hStripsTrkPlane"; string sPl=Form("%d",i); sName+=sPl; //cout<<"name="<<sName<<endl; TH1F* hStrip=new TH1F(sName.c_str(),sName.c_str(), 101,-1,100); hStrip->SetTitle(sName.c_str()); hStrip->GetXaxis()->SetTitle("Strip"); hStrip->GetXaxis()->CenterTitle(); hStrip->GetYaxis()->SetTitle(""); hStrip->GetYaxis()->CenterTitle(); hStrip->SetFillColor(0); //hStrip->SetBit(TH1::kCanRebin); vhStripsTrk[i]=hStrip; //strips are 0-95(96), pl 8,10: 4-67(64), pl 7,9: 0-63(64) //plane 6 has low strips in PI part (I think 0-67) //plane 11 has high strips in PI part (I think -95) //strips TPos: //plane 6 goes -2.64 -> 1.32 m //plane 11 goes -1.32 -> 2.64 m //plane 4,8,10 goes -2.40 -> 0.24 m //plane 5,7,9 goes -0.24 -> 2.40 m //2.64 - 2.40 = 24 cm = 5.85 strips //the FI planes "stick out" by ~6 strips sName="hTPosTrkPlane"; sName+=sPl; TH1F* hTPos=new TH1F(sName.c_str(),sName.c_str(),200,-4,4); hTPos->SetTitle(sName.c_str()); hTPos->GetXaxis()->SetTitle("TPos"); hTPos->GetXaxis()->CenterTitle(); hTPos->GetYaxis()->SetTitle(""); hTPos->GetYaxis()->CenterTitle(); hTPos->SetFillColor(0); //hTPos->SetBit(TH1::kCanRebin); vhTPosTrk[i]=hTPos; } //reconstruction object MeuReco reco; MeuSummary s; this->InitialiseLoopVariables(); //for(Int_t entry=0;entry<10000;entry++){ for(Int_t entry=0;entry<fEntries;entry++){ this->SetLoopVariables(entry); NtpStRecord& ntpst=(*fRec);//Make a reference instead of a pointer // this->FilterLI(ntpst,s); continue; TClonesArray& tcaTk=(*ntpst.trk); Int_t numTrks=tcaTk.GetEntries(); MAXMSG("MeuAnalysis",Msg::kInfo,10) <<"numTrks="<<numTrks <<", ndigit="<<fRec->evthdr.ndigit <<", nTrack="<<fRec->evthdr.ntrack<<endl; MeuSummary s; s.RFid=1.0; map<Int_t,MeuHitInfo> cp; hNTrack->Fill(numTrks); if (numTrks>1) continue; //Loop over tracks for (Int_t itrk=0;itrk<numTrks;itrk++){ const NtpSRTrack* ptrk= dynamic_cast<NtpSRTrack*>(tcaTk[itrk]); const NtpSRTrack& trk=(*ptrk); MAXMSG("MeuAnalysis",Msg::kInfo,10) <<"range="<<trk.range<<" nstrip="<<trk.nstrip<<" ds="<<trk.ds <<" fidend.dr,dz="<<trk.fidend.dr<<" "<<trk.fidend.dz<<" end-vtx.z:" <<trk.end.z-trk.vtx.z<<" "<<trk.stpz[1]-trk.stpz[0]<<endl; hStPerTrk->Fill(trk.nstrip); hYvsXVtx->Fill(trk.vtx.x,trk.vtx.y); hYvsXEnd->Fill(trk.end.x,trk.end.y); if (trk.end.plane>5 && trk.end.plane<115){ hYvsXVtxFid->Fill(trk.vtx.x,trk.vtx.y); hYvsXEndFid->Fill(trk.end.x,trk.end.y); } TClonesArray& tcaStp=(*ntpst.stp); Int_t numStps=tcaStp.GetEntries(); MAXMSG("MeuAnalysis",Msg::kDebug,10) <<"numStps="<<numStps<<endl; Float_t minR=999; Float_t maxR=0; for (Int_t i=0;i<trk.nstrip;i++){ //check for bug where strip index is -1 if (trk.stp[i]<0) { MAXMSG("MeuAnalysis",Msg::kInfo,500) <<"Skipping strip with trk.stp[i]="<<trk.stp[i]<<endl; continue; } //const NtpSRStrip* pstp=dynamic_cast<NtpSRStrip*>(tcaStp[i]); const NtpSRStrip* pstp= dynamic_cast<NtpSRStrip*>(tcaStp[trk.stp[i]]); const NtpSRStrip& stp=(*pstp); MAXMSG("MeuAnalysis",Msg::kInfo,20) <<"Strip="<<stp.strip<<", pl="<<stp.plane <<", sigCor="<<stp.ph1.sigcor<<", tpos="<<stp.tpos<<endl; vhStripsTrk[stp.plane]->Fill(stp.strip); vhTPosTrk[stp.plane]->Fill(stp.tpos); if (stp.plane>120) continue; hTPos->Fill(stp.tpos); hSigCor->Fill(stp.ph1.sigcor); hView->Fill(stp.planeview); hPlane->Fill(stp.plane); hStrip->Fill(stp.strip); hZ->Fill(stp.z); //fill the calibpos MeuHitInfo& c=cp[stp.plane]; c.Plane=stp.plane; c.View=stp.planeview; c.Strip=stp.strip; c.TPos=stp.tpos; c.SigCor2=stp.ph1.sigcor; c.X=trk.stpx[i]; c.Y=trk.stpy[i]; c.Z=trk.stpz[i]; hTPosVsPlane->Fill(stp.plane,stp.tpos); MAXMSG("MeuAnalysis",Msg::kInfo,10) <<"x="<<trk.stpx[i] <<", y="<<trk.stpy[i] <<", z="<<trk.stpz[i]<<endl; hYvsX->Fill(trk.stpx[i],trk.stpy[i]); hYvsZ->Fill(trk.stpz[i],trk.stpy[i]); hXvsZ->Fill(trk.stpz[i],trk.stpx[i]); hYvsPlane->Fill(stp.plane,trk.stpy[i]); hXvsPlane->Fill(stp.plane,trk.stpx[i]); Float_t r=sqrt(pow((trk.stpx[i]-1.5),2)+pow(trk.stpy[i],2)); if (r<minR) minR=r; if (r>maxR) maxR=r; hR->Fill(r); } //fill r histos if (maxR!=0 && minR!=999) { hMaxR->Fill(maxR); hMinR->Fill(minR); } } TClonesArray& tcaStp=(*ntpst.stp); Int_t numStps=tcaStp.GetEntries(); for (Int_t i=0;i<numStps;i++){ const NtpSRStrip* pstp=dynamic_cast<NtpSRStrip*>(tcaStp[i]); const NtpSRStrip& stp=(*pstp); //MAXMSG("MeuAnalysis",Msg::kInfo,200) // <<"Strip="<<stp.strip<<", pl="<<stp.plane<<endl; vhStrips[stp.plane]->Fill(stp.strip); vhTPos[stp.plane]->Fill(stp.tpos); } }//end of for MSG("MeuAnalysis",Msg::kInfo)<<"Finished main loop"<<endl; MSG("MeuAnalysis",Msg::kInfo) <<" ** Finished BasicPlots method **"<<endl; }

void MeuAnalysis::BasicReco (   )

<<PCCounter-bothSMCounter-coilHitCounter- Definition at line 2926 of file MeuAnalysis.cxx. References NtpSRPlane::beg, MeuReco::CalcFCPC(), MeuReco::CalcVtx(), MeuReco::CalcVtxSpecialND(), MeuReco::CalcWindow(), MeuReco::CheckWinContainment(), MeuSummary::DistToEdgeFid, NtpSRPlane::end, MeuSummary::EndPlane, MeuSummary::Event, NtpStRecord::evt, MeuCuts::ExtractPlInfo(), MeuSummary::FC, MeuCuts::FillSTSumDetails(), MeuCuts::FillSTSumRecoDetails(), RecRecordImp< T >::GetHeader(), MeuSummaryWriter::GetMeuSummaryToFill(), RecDataHeader::GetRun(), RecPhysicsHeader::GetSnarl(), InitialiseLoopVariables(), MAXMSG, MSG, NtpSRTrack::nstrip, OpenTxtFile(), MeuSummary::PC, NtpSRTrack::plane, Recalibrate(), MeuReco::Reconstruct(), MeuSummary::RFid, s(), SetLoopVariables(), MeuSummaryWriter::SummaryTreeFill(), MeuSummaryWriter::SummaryTreeFinish(), NtpStRecord::trk, MeuSummary::VtxPlane, and MeuSummary::WinSigMap. { MSG("MeuAnalysis",Msg::kInfo) <<" ** Running BasicReco method... **"<<endl; //open the output file for the histograms //fOutFile=this->OpenFile(100,"BasicReco"); MeuSummaryWriter summaryOut; TH1F* hWinSigMap=new TH1F("hWinSigMap","hWinSigMap",3000,-1,4000); hWinSigMap->SetTitle("hWinSigMap"); hWinSigMap->GetXaxis()->SetTitle("hWinSigMap"); hWinSigMap->GetXaxis()->CenterTitle(); hWinSigMap->SetFillColor(0); hWinSigMap->SetLineWidth(3); //hWinSigMap->SetBit(TH1::kCanRebin); Int_t FCCounter=0; Int_t PCCounter=0; Int_t TGCounter=0; Int_t badWindowCounter=0; Int_t badRecoCounter=0; Int_t goodWindowCounter=0; Int_t badFidCounter=0; ofstream& eventsTxtFile=*(this->OpenTxtFile(100,"myevents")); //reconstruction object MeuReco reco; MeuCuts cuts; //variable to turn on all the useful messages if required Msg::LogLevel_t logLevel=Msg::kDebug; this->InitialiseLoopVariables(); //for(Int_t entry=0;entry<10000;entry++){ for(Int_t entry=0;entry<fEntries;entry++){ this->SetLoopVariables(entry); MSG("MeuAnalysis",logLevel) <<"Entry="<<entry<<endl; //make a reference instead of a pointer NtpStRecord& ntp=(*fRec); const RecCandHeader& rec=ntp.GetHeader(); MAXMSG("MeuAnalysis",Msg::kInfo,1) <<"First entry: run="<<rec.GetRun() <<", snarl="<<rec.GetSnarl()<<endl; //MeuSummary s; MeuSummary& s=summaryOut.GetMeuSummaryToFill(); s.RFid=1.0; s.DistToEdgeFid=0.4; s.Event=entry; cuts.FillSTSumDetails(ntp,s,-1,-1); //a map to hold the info for each plane in the entry map<Int_t,MeuHitInfo> plInfo; //ensure there is only 1 track in the entry TClonesArray& trkTca=(*ntp.trk); Int_t numTrks=trkTca.GetEntries(); if (numTrks!=1) continue; const NtpSRTrack* ptrk=dynamic_cast<NtpSRTrack*>(trkTca[0]); const NtpSRTrack& trk=(*ptrk); //cut on number of planes if (sqrt(pow(1.*trk.plane.beg-trk.plane.end,2))<15) continue; MSG("MeuAnalysis",logLevel) <<"Entry="<<entry<<", run="<<rec.GetRun() <<", snarl="<<rec.GetSnarl() <<", trk.nstrips="<<trk.nstrip<<endl; TClonesArray& evtTca=(*ntp.evt); const NtpSREvent* pevt=dynamic_cast<NtpSREvent*>(evtTca[0]); cuts.ExtractPlInfo(ntp,s,plInfo,*pevt); reco.CalcVtx(s,plInfo); //check if the track extends into the calorimeter if (s.VtxPlane>120 && s.EndPlane>120) continue; //check if track stops in spectrometer if (s.EndPlane>120) continue; MSG("MeuAnalysis",logLevel)<<"CalcFCPC..."<<endl; reco.CalcFCPC(s,plInfo); //recalculate the vtx of the tracks coming in from the spectrometer MSG("MeuAnalysis",logLevel)<<"CalcVtxSpecialND..."<<endl; reco.CalcVtxSpecialND(s,plInfo); //cut on number of planes with potential new vertex if (sqrt(pow(1.*s.VtxPlane-s.EndPlane,2))<15) continue; Bool_t TG=!s.FC && !s.PC; if (s.PC) PCCounter++; else if (s.FC) FCCounter++; else if (TG) TGCounter++; if (s.FC && s.PC) MSG("MeuAnalysis",Msg::kError) <<"Both FC and PC!"<<endl; if (s.PC) { MSG("MeuAnalysis",logLevel)<<"Found PC event..."<<endl; Bool_t goodReco=reco.Reconstruct(s,plInfo); if (goodReco){ MSG("MeuAnalysis",logLevel)<<"Recalibrating..."<<endl; TClonesArray& evtTca=(*ntp.evt); const NtpSREvent* pevt=dynamic_cast<NtpSREvent*>(evtTca[0]); this->Recalibrate(ntp,s,plInfo,*pevt); MSG("MeuAnalysis",logLevel)<<"CalcWindow..."<<endl; reco.CalcWindow(s,plInfo); MSG("MeuAnalysis",logLevel)<<"CheckContainment..."<<endl; //fill the tree once event has been analysed if (s.WinSigMap<=0) { badWindowCounter++; } else if (!reco.CheckWinContainment(s,plInfo)){ badFidCounter++; //MSG("MeuAnalysis",Msg::kInfo)<<"Outside fid:"<<endl; } else { goodWindowCounter++; //Bool_t gapsOk=this->CheckTrackGaps(plInfo); //if (!gapsOk) { // MSG("MeuAnalysis",Msg::kInfo)<<"Large gap:"<<endl; // reco.PrintMeuHitInfo(plInfo); //} hWinSigMap->Fill(s.WinSigMap); eventsTxtFile<<entry<<endl; MAXMSG("MeuAnalysis",Msg::kInfo,10) <<"Good window: run=" <<rec.GetRun()<<", snarl="<<rec.GetSnarl() <<endl; cuts.FillSTSumRecoDetails(s,plInfo); summaryOut.SummaryTreeFill(plInfo); } } else { badRecoCounter++; MSG("MeuAnalysis",logLevel)<<"Bad reco..."<<endl; } } }//end of for MSG("MeuAnalysis",Msg::kInfo)<<"Finished main loop"<<endl; Double_t quantile=0.5;//quantile to compute Double_t meu=-1; hWinSigMap->GetQuantiles(1,&meu,&quantile); Double_t err=hWinSigMap->GetRMS()/sqrt(hWinSigMap->GetEntries()); Double_t relErr=err/meu; MSG("MeuAnalysis",Msg::kInfo) <<endl<<"*** Run Summary ***"<<endl <<"Raw MEU="<<meu<<" +/- "<<err<<" ("<<100.*relErr<<"%)"<<endl //<<endl <<"Total Entries="<<fEntries<<endl /*<<"Not LI="<<notLICounter<<endl <<"Not MM="<<notMMCounter<<endl <<"Good track="<<goodTrackCounter<<endl <<"Good length="<<goodLengthCounter<<endl <<"Good end & vtx="<<goodEndAndVtxCounter<<endl <<"Good XY reco="<<goodXYCounter<<endl <<endl <<"Good containmentPC="<<containmentCounter<<endl <<" both SM hit (when PC)="<<bothSMCounter <<" ("<<bothSM<<"% of PC)"<<endl <<" coil hit (when PC) ="<<coilHitCounter <<" ("<<coil_PC<<"% of PC)"<<endl <<" close LI (when PC) ="<<closeLICounter <<" ("<<closeLI<<"% of PC)"<<endl <<" dead Chips (when PC) ="<<deadChipsCounter <<" ("<<deadChips<<"% of PC)"<<endl <<" boundary TF (when PC)="<<boundaryTFCounter <<" ("<<boundaryTF<<"% of PC)"<<endl <<" bad views (when PC) ="<<badViewsCounter <<" ("<<badViews<<"% of PC)"<<endl <<" bad trk end (when PC)="<<badTrkEndCounter <<" ("<<badTrkEnd<<"% of PC)"<<endl <<" bad trk gap (when PC)="<<badTrkGapsCounter <<" ("<<badTrkGaps<<"% of PC)"<<endl <<" bad traceZ (when PC) ="<<badTraceZCounter <<" ("<<badTraceZ<<"% of PC)"<<endl <<" bad window (when PC) ="<<badWindowCounter <<" ("<<badWin<<"% of PC)"<<endl <<"Good containmentFC="<<containmentCounterFC<<endl <<" both SM hit (when FC)="<<bothSMCounterFC <<" ("<<bothSMFC<<"% of FC)"<<endl <<" coil hit (when FC) ="<<coilHitCounterFC <<" ("<<coil_FC<<"% of FC)"<<endl <<"Good containmentTG="<<containmentCounterTG<<endl <<" both SM hit (when TG)="<<bothSMCounterTG <<" ("<<bothSMTG<<"% of TG)"<<endl <<" coil hit (when TG) ="<<coilHitCounterTG <<" ("<<coil_TG<<"% of TG)"<<endl <<endl <<"Bad charge pos moment="<<momentCounter<<endl <<"Too high StripsPerPlane="<<stPerPlCounter<<endl <<" Either of above two="<<stPlOrMomentCounter<<endl <<endl <<"MC Info:"<<endl <<" Muon lowest energy >0.6 GeV="<<lowEnCounter<<endl <<" Muon lowest energy >5.0 GeV="<<lowEn5Counter<<endl <<endl */ <<"Containment summary:"<<endl <<" PC="<<PCCounter<<endl // <<" (PC - all cuts=" //closeLICounter-deadChipsCounter-boundaryTFCounter- //badViewsCounter-badTrkEndCounter-badTrkGapsCounter-badTraceZCounter- //badWindowCounter<<")"<<endl <<" FC="<<FCCounter<<endl//<<" (FC-coil-bothSM=" //<<FCCounter-coilHitCounterFC-bothSMCounterFC<<")"<<endl <<" TG="<<TGCounter<<endl//<<" (TG-coil-bothSM=" // <<TGCounter-coilHitCounterTG-bothSMCounterTG<<")"<<endl <<" FC+PC+TG="<<FCCounter+PCCounter+TGCounter<<endl <<endl <<"Total PC with good track window = "<<goodWindowCounter<<endl <<" bad track window = "<<badWindowCounter<<endl <<" bad fiducial vol = "<<badFidCounter<<endl <<" bad reco = "<<badRecoCounter<<endl <<endl; //finish the summary tree and output it summaryOut.SummaryTreeFinish(); if (fRec) delete fRec; MSG("MeuAnalysis",Msg::kInfo) <<" ** Finished BasicReco method **"<<endl; }

void MeuAnalysis::CheckTrackGaps2 (  const NtpStRecord &  ntp  )  const [private]

Definition at line 2087 of file MeuAnalysis.cxx. References MAXMSG, NtpSRTrack::nstrip, NtpSRStrip::plane, NtpSRTrack::stp, NtpStRecord::stp, and NtpStRecord::trk. { TClonesArray& tcaTk=(*ntp.trk); Int_t numTrks=tcaTk.GetEntries(); //Loop over tracks for (Int_t itrk=0;itrk<numTrks;itrk++){ const NtpSRTrack* ptrk= dynamic_cast<NtpSRTrack*>(tcaTk[itrk]); const NtpSRTrack& trk=(*ptrk); TClonesArray& tcaStp=(*ntp.stp); //this loop is just over the tracked strips for (Int_t i=0;i<trk.nstrip;i++){ //check for bug where strip index is -1 if (trk.stp[i]<0) { MAXMSG("MeuAnalysis",Msg::kInfo,500) <<"Skipping strip with trk.stp[i]="<<trk.stp[i]<<endl; continue; } const NtpSRStrip* pstp= dynamic_cast<NtpSRStrip*>(tcaStp[trk.stp[i]]); const NtpSRStrip& stp=(*pstp); //MAXMSG("MeuAnalysis",Msg::kInfo,200) // <<"Strip="<<stp.strip<<", pl="<<stp.plane // <<", sigCor="<<stp.ph1.sigcor<<", tpos="<<stp.tpos<<endl; Int_t pl=stp.plane; if (pl>120) continue; } } }

Bool_t MeuAnalysis::FilterLI (  const AtmosEvent &  ev  )  const [private]

Definition at line 4170 of file MeuAnalysis.cxx. References MAXMSG, AtmosStrip::Plane, AtmosStrip::QPEcorr, and AtmosEvent::StripList. { Float_t sumSigCorEast=0; Float_t sumSigCorWest=0; map<Int_t,Int_t> hpp; TClonesArray& strips=(*ev.StripList); Int_t numStrips=strips.GetEntries(); if (numStrips<=0) return false;//pass the event, not LI //loop over the strips in the snarl for (Int_t hit=0;hit<numStrips;hit++){ const AtmosStrip* strip=dynamic_cast<AtmosStrip*>(strips[hit]); const AtmosStrip& st=(*strip); sumSigCorEast+=st.QPEcorr[0]; sumSigCorWest+=st.QPEcorr[1]; hpp[st.Plane]++; }//end of loop over strips Float_t avHpp=0; vector<Float_t> pb(8);//to store fractions hit in each pb region //loop over all the planes flashed for (map<Int_t,Int_t>::iterator hppIt=hpp.begin(); hppIt!=hpp.end();++hppIt){ avHpp+=hppIt->second; Int_t pl=hppIt->first; //work out the fraction of the planes flashed by each pb //are actually hit if (pl>=1 && pl<=64) pb[0]+=1./64;//0=bookend else if (pl>=65 && pl<=128) pb[1]+=1./64; else if (pl>=129 && pl<=192) pb[2]+=1./64; else if (pl>=193 && pl<=248) pb[3]+=1./56;//249=bookend else if (pl>=250 && pl<=313) pb[4]+=1./64; else if (pl>=314 && pl<=377) pb[5]+=1./64; else if (pl>=378 && pl<=441) pb[6]+=1./64; else if (pl>=442 && pl<=485) pb[7]+=1./44; } if (hpp.size()>0) avHpp/=hpp.size(); Float_t sumSig=sumSigCorEast+sumSigCorWest; Float_t minSigCorEW = TMath::Min(sumSigCorEast, sumSigCorWest); Float_t asym=0; if (sumSig) asym = (sumSig - 2 * minSigCorEW) / sumSig; Float_t fractFlashed1=0; Float_t fractFlashed2=0; for (vector<Float_t>::iterator it=pb.begin(); it!=pb.end();++it){ Float_t fract=*it; if (fract>fractFlashed1){ //copy the 1 into 2 fractFlashed2=fractFlashed1; //then set new highest fractFlashed fractFlashed1=fract; } MAXMSG("MeuAnalysis",Msg::kVerbose,200) <<"fract="<<fract<<", f1="<<fractFlashed1 <<", f2="<<fractFlashed2<<endl; } Float_t ratioFlashed=0; if (fractFlashed1>0) ratioFlashed=fractFlashed2/fractFlashed1; Bool_t li=false; //check if LI if ( (asym>0.5 || minSigCorEW>1.7e6) && avHpp>3 && ratioFlashed<0.05 && fractFlashed1>0.85 && fractFlashed2<0.1 ){ li=true; MAXMSG("MeuAnalysis",Msg::kDebug,200) <<"LI Event:"<<endl <<" Asymmetry="<<asym<<endl <<" Av. HPP ="<<avHpp<<endl <<" Fract Fl1="<<fractFlashed1<<", Fract Fl2="<<fractFlashed2 <<" Ratio Fla="<<ratioFlashed<<endl; } return li; }

Bool_t MeuAnalysis::FilterLI (  const NtpStRecord &  ntp, const MeuSummary &  s )  const [private]

Definition at line 4025 of file MeuAnalysis.cxx. References MeuSummary::Detector, MAXMSG, NtpSRStrip::ph0, NtpSRStrip::ph1, NtpSRStrip::plane, s(), NtpSRPulseHeight::sigcor, NtpStRecord::stp, NtpSRStrip::time0, and NtpSRStrip::time1. Referenced by BasicPlots(), and MakeSummaryTreeWithAtNu(). { //only analyse if FD if (s.Detector!=Detector::kFar) return false; const Bool_t cleanTheEvent=false; Float_t sumSigCorEast=0; Float_t sumSigCorWest=0; map<Int_t,Int_t> hpp; // <plane, # of hits> multiset<Double_t> times; Int_t numStrips=ntp.stp->GetEntries(); if (numStrips<=0) return false;//pass the event, not LI //loop over the strips in the snarl for (Int_t hit=0;hit<numStrips;hit++){ NtpSRStrip* strip = dynamic_cast<NtpSRStrip*>((*ntp.stp)[hit]); Double_t time=strip->time0; if (time<0 || time>1) time=strip->time1; if (time>0 && time<=1) { times.insert(time); } //else just don't put the time in the map - clearly rubbish }//end of loop over strips //get the median time from the map multiset<Double_t>::const_iterator it=times.begin(); advance(it,times.size()/2); Double_t medianTime=*it; MAXMSG("MeuAnalysis",Msg::kDebug,100) <<"FilterLI: Median time="<<medianTime<<endl; //second loop for (Int_t hit=0;hit<numStrips;hit++){ NtpSRStrip* strip = dynamic_cast<NtpSRStrip*>((*ntp.stp)[hit]); if (cleanTheEvent) { if (strip->ph0.sigcor > 0){ if (strip->time0<medianTime-200e-9 || strip->time0>medianTime+200e-9) continue; } if (strip->ph1.sigcor > 0){ if (strip->time1<medianTime-200e-9 || strip->time1>medianTime+200e-9) continue; } } // if (cleanTheEvent) sumSigCorEast+=strip->ph0.sigcor; sumSigCorWest+=strip->ph1.sigcor; hpp[strip->plane] += 1; }//end of loop over strips Float_t avHpp=0; vector<Float_t> pb(8);//to store fractions hit in each pb region //loop over all the planes flashed for (map<Int_t,Int_t>::iterator hppIt=hpp.begin(); hppIt!=hpp.end();++hppIt){ avHpp+=hppIt->second; Int_t pl=hppIt->first; if (pl>=1 && pl<=64) pb[0]+=1./64;//0=bookend else if (pl>=65 && pl<=128) pb[1]+=1./64; else if (pl>=129 && pl<=192) pb[2]+=1./64; else if (pl>=193 && pl<=248) pb[3]+=1./56;//249=bookend else if (pl>=250 && pl<=313) pb[4]+=1./64; else if (pl>=314 && pl<=377) pb[5]+=1./64; else if (pl>=378 && pl<=441) pb[6]+=1./64; else if (pl>=442 && pl<=485) pb[7]+=1./44; MAXMSG("MeuAnalysis",Msg::kDebug,2000) <<"plane="<<pl<<", hpp="<<hppIt->second<<endl; } if (hpp.size()>0) avHpp/=hpp.size(); Float_t sumSig=sumSigCorEast+sumSigCorWest; Float_t minSigCorEW = TMath::Min(sumSigCorEast, sumSigCorWest); Float_t asym=0; if (sumSig) asym = (sumSig - 2 * minSigCorEW) / sumSig; Float_t fractFlashed1=0; Float_t fractFlashed2=0; for (vector<Float_t>::iterator it=pb.begin(); it!=pb.end();++it){ Float_t fract=*it; if (fract>fractFlashed1){ //copy the 1 into 2 fractFlashed2=fractFlashed1; //then set new highest fractFlashed fractFlashed1=fract; } //make sure to get fractFlashed2 if fractFlashed2<fract<frFla1 else if (fract>fractFlashed2) { fractFlashed2=fract; } MAXMSG("MeuAnalysis",Msg::kDebug,200) <<"fract="<<fract<<", f1="<<fractFlashed1 <<", f2="<<fractFlashed2<<endl; } Float_t ratioFlashed=0; if (fractFlashed1>0) ratioFlashed=fractFlashed2/fractFlashed1; Bool_t li=false; // check if LI // don't check the asymmetry if the VA is saturating if ( (asym>0.5 || minSigCorEW>1.7e6) && avHpp>3 && ratioFlashed<0.05 && fractFlashed1>0.85 && fractFlashed2<0.1 ){ li=true; MAXMSG("MeuAnalysis",Msg::kInfo,10) <<"LI Event:"<<endl <<" Asymmetry="<<asym<<endl <<" Av. HPP ="<<avHpp<<endl <<" Fract Fl1="<<fractFlashed1 <<", Fract Fl2="<<fractFlashed2 <<" Ratio Fla="<<ratioFlashed<<endl; } else{ MAXMSG("MeuAnalysis",Msg::kDebug,100) <<"Non-LI Event:"<<endl <<" Asymmetry="<<asym<<endl <<" Av. HPP ="<<avHpp<<endl <<" Fract Fl1="<<fractFlashed1 <<", Fract Fl2="<<fractFlashed2 <<" Ratio Fla="<<ratioFlashed<<endl; } return li; }

std::string MeuAnalysis::GetFirstFileName (  std::string  wildcardString  )  const

Definition at line 375 of file MeuAnalysis.cxx. References Form(), gSystem(), and s(). Referenced by MakeChain(). { //check if there is actually a wildcard if (!TString(wildcardString.c_str()).MaybeWildcard()) { return wildcardString; } //this was ripped off from ROOT's TChain::AddFile vector<string> fileList; //wildcarding used in name TString basename(wildcardString.c_str()); Int_t dotslashpos = basename.Index(".root/"); TString behind_dot_root; if (dotslashpos>=0) { // Copy the tree name specification behind_dot_root = basename(dotslashpos+6,basename.Length()-dotslashpos+6); // and remove it from basename basename.Remove(dotslashpos+5); } Int_t slashpos = basename.Last('/'); TString directory; if (slashpos>=0) { directory = basename(0,slashpos); // Copy the directory name basename.Remove(0,slashpos+1); // and remove it from basename } else { directory = gSystem->WorkingDirectory(); } const char *file; void *dir = gSystem->OpenDirectory(gSystem->ExpandPathName(directory.Data())); if (dir) { //create a TList to store the file names (not yet sorted) TList l; TRegexp re(basename,kTRUE); while ((file = gSystem->GetDirEntry(dir))) { if (!strcmp(file,".") || !strcmp(file,"..")) continue; TString s = file; if ( (basename!=file) && s.Index(re) == kNPOS) continue; l.Add(new TObjString(file)); } gSystem->FreeDirectory(dir); //sort the files in alphanumeric order l.Sort(); TIter next(&l); TObjString *obj; while ((obj = (TObjString*)next())) { file = obj->GetName(); if (behind_dot_root.Length()!=0){ string fileName=Form("%s/%s/%s",directory.Data(), file,behind_dot_root.Data()); fileList.push_back(fileName); //nf += AddFile(Form("%s/%s/%s",directory.Data(),file,behind_dot_root.Data()),kBigNumber); } else { string fileName=Form("%s/%s",directory.Data(),file); fileList.push_back(fileName); //nf += AddFile(Form("%s/%s",directory.Data(),file),kBigNumber); } } l.Delete(); } //check if any files were found if (fileList.begin()!=fileList.end()){ cout<<"Used wildcard expansion to find first file name=" <<*fileList.begin()<<endl; return *fileList.begin(); } else { return "";//null string } }

Float_t MeuAnalysis::GetTemperature (  const MeuSummary &  s  )  const [private]

Definition at line 4261 of file MeuAnalysis.cxx. References MeuSummary::Detector, MAXMSG, MSG, VldTimeStamp::Print(), s(), and MeuSummary::TimeSec. Referenced by MakeSummaryTreeWithAtNu(). { static Bool_t firstCall=true; static Bool_t foundFile=false; static vector <Float_t> vTemperatures; static Int_t firstTime=-1; static Int_t lastTime=-1; //only get the temperature for the FD at present if (s.Detector!=Detector::kFar) { MAXMSG("MeuAnalysis",Msg::kInfo,1) <<"Not retrieving the temperature for dt="<<s.Detector<<endl; return 0; } if (firstCall){ char* envVariable=getenv("MEUDATA"); if (envVariable==NULL){ MAXMSG("MeuAnalysis",Msg::kWarning,1) <<"Can't find ENV variable \"MEUDATA\" to" <<" open temperatures file"<<endl; firstCall=false; return 0; } string sEnv=envVariable; string sTemperatureFileName=sEnv+"/temps.dat"; MSG("MeuAnalysis",Msg::kInfo) <<"GetTemperature: Looking for file: " <<sTemperatureFileName<<endl; ifstream temperatureFile(sTemperatureFileName.c_str()); if(!temperatureFile){ MAXMSG("MeuAnalysis",Msg::kWarning,1) <<"Can't find file of temperatures"<<endl; firstCall=false; foundFile=false; return 0; } foundFile=true; map<Int_t,Float_t> temperatures; MSG("MeuAnalysis",Msg::kInfo) <<"Loading temperatures from txt file..."<<endl; if (temperatureFile){ Float_t Tf=-1; Float_t Tc=-1; Int_t time=0; while(temperatureFile>>time>>Tf) { Tc=(Tf-32)/1.8; temperatures[time]=Tc; } } firstTime=temperatures.begin()->first; lastTime=(--temperatures.end())->first; MSG("MeuAnalysis",Msg::kInfo) <<"Found first time="<<firstTime<<", lastTime="<<lastTime<<endl; VldTimeStamp first(firstTime,0); first.Print(); VldTimeStamp last(lastTime,0); last.Print(); Int_t lastIndex=(lastTime-firstTime)/300; MSG("MeuAnalysis",Msg::kInfo)<<"Size of table="<<lastIndex<<endl; vTemperatures.reserve(lastIndex+1); for (Int_t i=0;i<lastIndex+1;i++) vTemperatures.push_back(0); //loop and fill vector for (map<Int_t,Float_t>::iterator it=temperatures.begin(); it!=temperatures.end();++it){ Int_t time=it->first; Int_t index=(time-firstTime)/300;//every 5 mins MAXMSG("MeuAnalysis",Msg::kDebug,50) <<"time="<<it->first <<", mins since first="<<(time-firstTime)/300 <<", temperature="<<it->second<<endl; vTemperatures[index]=it->second; } } firstCall=false; if (!foundFile) return 0; Int_t ind=(s.TimeSec-firstTime)/300; Float_t Tc=0; Int_t lastIndex=(lastTime-firstTime)/300; Int_t timeSkipped=0; if (s.TimeSec>=firstTime && s.TimeSec<=lastTime){ while (ind<=lastIndex && ind>0){ MAXMSG("MeuAnalysis",Msg::kVerbose,500) <<"500 mins T="<<vTemperatures[ind]<<", ind="<<ind<<endl; if (vTemperatures[ind]>0){ MAXMSG("MeuAnalysis",Msg::kDebug,50) <<"Found good T="<<vTemperatures[ind]<<endl; Tc=vTemperatures[ind]; break; } timeSkipped++; ind++; } } else MSG("MeuAnalysis",Msg::kWarning)<<"Time is out of range"<<endl; //check that not too big a time (1 hour) was skipped if (timeSkipped<=12) return Tc; else { MAXMSG("MeuAnalysis",Msg::kInfo,50) <<"Too big a gap (>1 hour) in temperature reading" <<" to be reliable ("<<timeSkipped<<"*5 mins)." <<" Returning 0 degC"<<endl; return 0; } }

void MeuAnalysis::InitialiseLoopVariables (   )  [private]

Definition at line 762 of file MeuAnalysis.cxx. References MSG. Referenced by BasicPlots(), BasicReco(), MakeSummaryTreeWithAtNu(), MakeSummaryTreeWithNtpSt(), N_1Plots(), SnarlList(), SpillPlots(), and TestEventLoop(). { MSG("MeuAnalysis",Msg::kDebug)<<"Initialising loop variables..."<<endl; MSG("MeuAnalysis",Msg::kDebug)<<"Initialisation complete"<<endl; }

void MeuAnalysis::InputFileName (  std::string  f  )  [static]

Definition at line 153 of file MeuAnalysis.cxx. References fInputFileName, and MSG. { if (f!=""){ MSG("MeuAnalysis",Msg::kInfo) <<"Running with input file name="<<f<<endl; fInputFileName=f; } }

void MeuAnalysis::LoadTrees (  Bool_t  loadTrees  )  [static]

Definition at line 164 of file MeuAnalysis.cxx. References fLoadTrees, and MSG. { MSG("MeuAnalysis",Msg::kInfo) <<"Setting fLoadTrees="<<loadTrees<<endl; fLoadTrees=loadTrees; }

void MeuAnalysis::MakeChain (   )  [private]

Definition at line 455 of file MeuAnalysis.cxx. References fChain, fChainBD, fUseAtNu, GetFirstFileName(), MakeFileList(), MSG, and ObjectExistsInFile(). { //get the files to open vector<string> fileList=this->MakeFileList(); TFile* firstFile=0; // case with one single file if (!TString(*fileList.begin()).MaybeWildcard()) { firstFile=TFile::Open((*fileList.begin()).c_str(),"READ"); } else { MSG("MeuAnalysis",Msg::kInfo) <<"Found a wildcard present in list of file names: " <<*fileList.begin() <<endl<<"Looking for first filename using wildcard..."<<endl; string fileName=this->GetFirstFileName(*fileList.begin()); if (fileName!="") { firstFile=TFile::Open(fileName.c_str(),"READ"); } else { MSG("MeuAnalysis",Msg::kInfo) <<"No files found matching:"<<endl <<*fileList.begin()<<", will exit here..."<<endl; exit(1); } } //create a chain if (false && fUseAtNu) { if (this->ObjectExistsInFile(firstFile,"ntp")) { fChain=new TChain("ntp"); } } else { //check if the AtmosEvent ntuple exists in the file if (this->ObjectExistsInFile(firstFile,"ntp")) { fChain=new TChain("ntp"); fUseAtNu=true;//set this to true } //else check if the SR ntuple exists else if (this->ObjectExistsInFile(firstFile,"NtpSt")) { fChain=new TChain("NtpSt"); } else { MSG("MeuAnalysis",Msg::kInfo) <<"Not creating NtpSt TChain because it does not exist in file" <<endl; } //check in BD ntuple exists if (this->ObjectExistsInFile(firstFile,"NtpBDLite")) { fChainBD=new TChain("NtpBDLite");//have to load libraries too } else { MSG("MeuAnalysis",Msg::kInfo) <<"Not creating NtpBDLite TChain because it does not exist" <<" in file"<<endl; } } firstFile->Close(); if (firstFile) delete firstFile; firstFile=0; Int_t nf=0; Int_t nfBD=0; //add the files to the chain for (vector<string>::iterator file=fileList.begin(); file!=fileList.end();++file){ //test if file already exists ifstream openOk((*file).c_str()); //check if a wildcard was used because ifstream can't open wildcards Int_t stars=(*file).find("*"); Int_t quest=(*file).find("?"); //check if file existed if (!openOk && !(quest>=0 || stars>=0)){ MSG("MeuAnalysis",Msg::kInfo) <<endl<<endl <<"***********************************************************" <<endl<<"Can't find file="<<*file<<endl <<"Note: you can't use '~/'. It has to be the full path"<<endl <<"***********************************************************" <<endl<<endl <<"Exiting here!"<<endl; exit(0); } MSG("MeuAnalysis",Msg::kInfo)<<"Adding file(s)="<<*file<<endl; nf+=fChain->Add((*file).c_str()); if (fChainBD) nfBD+=fChainBD->Add((*file).c_str()); //fChain->Print();//very verbose } if(nf==0){ MSG("MeuAnalysis",Msg::kFatal) <<endl<<endl <<"*************************************************************" <<endl<<"No *.root files found"<<endl <<"Please set MEUDATA to the directory containing the" <<" *.root files"<<endl <<"Or give the txt file containing the files to be input"<<endl <<"Note: If more than one file is found they will be" <<" concatenated in a TChain and treated as one"<<endl <<"*************************************************************" <<endl<<endl<<"Program will exit here"<<endl; exit(0); } MSG("MeuAnalysis",Msg::kInfo) <<endl<<"Added "<<nf<<" file(s) to TChain"<<endl; MSG("MeuAnalysis",Msg::kInfo) <<endl<<"Added "<<nfBD<<" beam data (BD) file(s) to TChain" <<endl; if (fChain) { MSG("MeuAnalysis",Msg::kInfo) <<"Ntuple information:"<<endl; fChain->Show(0); if (fChainBD) { MSG("MeuAnalysis",Msg::kInfo) <<"NtpBDLite information:"<<endl; fChainBD->Show(0); } } MSG("MeuAnalysis",Msg::kInfo) <<endl<<"Added "<<nfBD<<" beam data (BD) file(s) to TChain." <<endl; MSG("MeuAnalysis",Msg::kInfo) <<endl<<"Analysing "<<nf<<" file(s). Reading from disk..."<<endl; }

vector< string > MeuAnalysis::MakeFileList (   )  [private]

Check the fInputFileName first then check the env variable Definition at line 259 of file MeuAnalysis.cxx. References fInputFileName, and MSG. Referenced by MakeChain(). { vector<string> fileList; if (fInputFileName!=""){ Int_t findGives=fInputFileName.find(".root"); MSG("MeuAnalysis",Msg::kInfo) <<"Checking input string for \"*.root\", position in string=" <<findGives<<endl; if (findGives>0){ //add the file direct to the list fileList.push_back(fInputFileName); } else{//treat the file as a list of root files ifstream inputFile(fInputFileName.c_str()); //check if file exists if (inputFile){ //variables to hold input from file string file=""; //read in from the text file and fill objects while(inputFile>>file) { MSG("MeuAnalysis",Msg::kDebug) <<"Found input file name="<<file<<endl; fileList.push_back(file); } MSG("MeuAnalysis",Msg::kDebug) <<"Files names found in txt file="<<fileList.size()<<endl; } else{ MSG("MeuAnalysis",Msg::kFatal) <<endl<<endl <<"**********************************************************" <<endl<<"Input txt file of file names does not exist!"<<endl <<"InputFileName="<<fInputFileName<<endl <<"**********************************************************" <<endl<<endl<<"Program will exit here"<<endl; exit(0); } } } //return the fileList if files were found if (fileList.size()>0) return fileList; //Check the env variable to find files char* envVariable=getenv("MEUDATA"); if (envVariable==NULL){ MSG("MeuAnalysis",Msg::kFatal) <<endl<<endl <<"*************************************************************" <<endl<<"Environmental variable MEUDATA not set!"<<endl <<"Please set MEUDATA to the directory containing the" <<" ntuple root files"<<endl <<"Note: If more than one file is found they will be" <<" concatenated and treated as one"<<endl <<"*************************************************************" <<endl<<endl<<"Program will exit here"<<endl; exit(0); } string sEnv=""; sEnv=envVariable; sEnv+="/*.root"; MSG("MeuAnalysis",Msg::kInfo) <<"Looking for *.root files using the env variable"<<endl <<"MEUDATA="<<sEnv<<endl; fileList.push_back(sEnv); return fileList; }

void MeuAnalysis::MakeSummaryTree (   )

Definition at line 4380 of file MeuAnalysis.cxx. References MakeSummaryTreeWithAtNu(), and MakeSummaryTreeWithNtpSt(). { //this method just adds a level of indirection so that you only //need a single macro to make the summary tree if (fUseAtNu){ this->MakeSummaryTreeWithAtNu(); } else{ this->MakeSummaryTreeWithNtpSt(); } }

void MeuAnalysis::MakeSummaryTreeWithAtNu (   )  [private]

Definition at line 4395 of file MeuAnalysis.cxx. References MeuCuts::AnalyseCoilProximity(), AtmosTrack::AtNuNplanes, MeuCounter::bad2ndContCounter, MeuCounter::badDistEndCounter, MeuCounter::badEndGapsCounter, MeuCounter::badFidCounter, MeuCounter::badRecoCounter, MeuCounter::badTraceZCounter, MeuCounter::badTrackEndCounter, MeuCounter::badViewsCounter, MeuCounter::badWindowCounter, MeuCounter::bothSMHitCounter, MeuReco::CalcFCPC(), MeuReco::CalcStripDists(), MeuReco::CalcTrace(), MeuReco::CalcVtx(), MeuReco::CalcWindow(), MeuCuts::CheckTrkViews(), MeuReco::CheckWinContainment(), MeuCounter::coilHitCounter, MeuSummary::DistToEdgeFid, MeuSummary::EndDistToEdge, MeuSummary::EndPlane, AtmosTrack::EndTraceZ, MeuSummary::Event, MeuCuts::ExtractPlInfo(), MeuCuts::ExtractTruthInfo(), MeuSummary::FC, MeuCounter::FCCounter, MeuCuts::FillEnVsDist(), MeuCuts::FillEnVsDist2(), MeuHistos::FillMeuHistos(), MeuCuts::FillSTSumDetails(), MeuCuts::FillSTSumRecoDetails(), MeuCuts::FilterBadDistEndStrip(), MeuCuts::FilterBadEndGaps(), MeuCuts::FilterBadTrackEnd(), MeuCuts::FilterBadXY(), FilterLI(), fRun, fSubrun, MeuSummaryWriter::GetMeuSummaryToFill(), GetTemperature(), MeuCounter::goodContainmentCounter, MeuCounter::goodWindowCounter, MeuCounter::goodXYCounter, InitialiseLoopVariables(), MeuCounter::longTrackCounter, MeuHitInfo::LPos, MAXMSG, MeuSummary::MCHighEn, MeuSummary::MCLowEn, MSG, MeuCounter::notLICounter, AtmosEvent::NTracks, MeuCounter::oneTrackCounter, MeuSummary::PC, MeuCounter::PCCounter, MeuHitInfo::Plane, print(), MeuHistos::PrintMeuValues(), MeuCounter::PrintNtpSt(), RecalibrateAtNu(), MeuReco::Reconstruct(), MeuSummary::RFid, MeuSummary::RFidCoil, MeuSummary::Run, AtmosEvent::Run, s(), SetLoopVariables(), MeuSummary::SMBoth, MeuSummary::Snarl, AtmosEvent::Snarl, MeuSummaryWriter::SummaryTreeFill(), MeuSummaryWriter::SummaryTreeFinish(), MeuSummaryWriter::SummaryTreeSetup(), MeuSummary::Temperature, MeuCounter::TGCounter, MeuHitInfo::TPos, AtmosEvent::TrackList, MeuSummary::VtxDistToEdge, MeuSummary::VtxPlane, MeuSummary::WinSigMap, MeuHitInfo::X, MeuHitInfo::Y, and MeuHitInfo::Z. Referenced by MakeSummaryTree(). { MSG("MeuAnalysis",Msg::kInfo) <<" ** Running MakeSummaryTreeWithAtNu method... **"<<endl; //open the output file for the histograms //fOutFile=this->OpenFile(100,"MakeSummaryTreeFar"); MeuSummaryWriter summaryOut; summaryOut.SummaryTreeSetup(fRun,fSubrun); //string temps="myeventsLowEn"; //if (fRun>10000) temps+="MC"; //ofstream& eventsTxtFile=*(this->OpenTxtFile(100,"myevents")); //ofstream& txtFileLowEn=*(this->OpenTxtFile(100,temps.c_str())); //ofstream& txtFileNonPC=*(this->OpenTxtFile(100,"myeventsNonPC")); //reconstruction object const MeuReco reco; const MeuCuts cuts; const MeuHistos histos; MeuCounter cnt; //variable to turn on all the useful messages if required Msg::LogLevel_t logLevel=Msg::kDebug; this->InitialiseLoopVariables(); //for(Int_t entry=24196;entry<fEntries;entry++){ //for(Int_t entry=20307;entry<fEntries;entry++){//crashed on 32998 //for(Int_t entry=30000;entry<fEntries;entry++){//didn't crash 32998 //for(Int_t entry=0;entry<2000;entry++){ for(Int_t entry=0;entry<fEntries;entry++){ this->SetLoopVariables(entry); MSG("MeuAnalysis",logLevel) <<"Entry="<<entry<<endl; //make a reference instead of a pointer AtmosEvent& ev=(*fev); MAXMSG("MeuAnalysis",Msg::kInfo,1) <<"First entry: run="<<ev.Run <<", snarl="<<ev.Snarl<<endl; //MSG("MeuAnalysis",Msg::kInfo) //<<"run="<<ev.Run //<<", snarl="<<ev.Snarl //<<", entry="<<entry //<<", unixTime="<<ev.UnixTime //<<", #scintHits="<<ev.NScintHits //<<endl; //cut out the LI events if (this->FilterLI(ev)){ //eventsTxtFileLI<<event<<endl; continue; } cnt.notLICounter++; //cut out events with multiple tracks if (ev.NTracks!=1) continue; cnt.oneTrackCounter++; const AtmosTrack* track=dynamic_cast<const AtmosTrack*> (ev.TrackList->At(0)); const AtmosTrack& trk=(*track); //cut out events with less than 20 planes Bool_t goodLength=trk.AtNuNplanes>20; if (!goodLength){ continue; } cnt.longTrackCounter++; //MeuSummary s; MeuSummary& s=summaryOut.GetMeuSummaryToFill(); s.RFid=3.0*Munits::m; s.RFidCoil=0.5*Munits::m; s.DistToEdgeFid=1.0*Munits::m; s.Event=entry; cuts.FillSTSumDetails(ev,s); //a map to hold the info for each plane in the entry map<Int_t,MeuHitInfo> plInfo; cuts.ExtractPlInfo(ev,s,plInfo); reco.CalcVtx(s,plInfo); MSG("MeuAnalysis",logLevel)<<"CalcFCPC..."<<endl; reco.CalcFCPC(s,plInfo); //check that the XY position is sensible if (cuts.FilterBadXY(ev,s)){ continue; } cnt.goodXYCounter++; Bool_t TG=!s.FC && !s.PC; if (s.PC) cnt.PCCounter++; else if (s.FC) cnt.FCCounter++; else if (TG) cnt.TGCounter++; if (s.FC && s.PC) MSG("MeuAnalysis",Msg::kError) <<"Both FC and PC!"<<endl; Bool_t awayFromCoil=cuts.AnalyseCoilProximity(ev,s); if (s.PC && !awayFromCoil) cnt.coilHitCounter++; if (s.PC && s.SMBoth && awayFromCoil) cnt.bothSMHitCounter++; Bool_t goodContainment=s.PC && awayFromCoil && !s.SMBoth; if (!goodContainment) continue; MSG("MeuAnalysis",logLevel)<<"Found PC event..."<<endl; cnt.goodContainmentCounter++; Int_t diffAtEnd=-1; Int_t diffAtVtx=-1; if (!cuts.CheckTrkViews(s,&diffAtVtx,&diffAtEnd)) { cnt.badViewsCounter++; continue; } if (cuts.FilterBadTrackEnd(s,plInfo)){ cnt.badTrackEndCounter++; continue; } Bool_t goodReco=reco.Reconstruct(s,plInfo); if (!goodReco){ MSG("MeuAnalysis",logLevel)<<"Bad reco..."<<endl; cnt.badRecoCounter++; continue; } if (cuts.FilterBadEndGaps(s,plInfo)){//after reco only cnt.badEndGapsCounter++; continue; } MSG("MeuAnalysis",logLevel)<<"Recalibrating..."<<endl; this->RecalibrateAtNu(ev,s,plInfo); MSG("MeuAnalysis",logLevel)<<"CalcWindow..."<<endl; reco.CalcWindow(s,plInfo); //MSG("MeuAnalysis",logLevel)<<"CheckContainment..."<<endl; //fill the tree once event has been analysed if (s.WinSigMap<=0) { cnt.badWindowCounter++; continue; } if (!reco.CheckWinContainment(s,plInfo)){//nothing for FD cnt.badFidCounter++; continue; } //recalculate the containment! //s.DistToEdgeFid=0.35;//allow a slight shift in position //reco.CalcFCPC(s,plInfo); if (!s.PC){ cnt.bad2ndContCounter++; static Bool_t print=false; if (print){ MAXMSG("MeuAnalysis",Msg::kInfo,100) <<endl<<endl<<"Recalculated containment and not PC!" <<" PC="<<s.PC<<", FC="<<s.FC <<" entry="<<entry <<", snarl="<<s.Snarl<<", run="<<s.Run<<endl <<"Low En="<<s.MCLowEn<<" high="<<s.MCHighEn<<endl <<endl; MAXMSG("MeuAnalysis",Msg::kInfo,100) <<"vtx Pl="<<s.VtxPlane<<", end Pl="<<s.EndPlane <<", dVtx="<<s.VtxDistToEdge <<", dEnd="<<s.EndDistToEdge<<endl; for (map<Int_t,MeuHitInfo>::const_iterator it= plInfo.begin(); it!=plInfo.end();++it){ const MeuHitInfo& c=it->second; MAXMSG("MeuAnalysis",Msg::kInfo,500) <<"pl="<<c.Plane<<", z="<<c.Z <<", x="<<c.X<<", y="<<c.Y <<" t="<<c.TPos<<", l="<<c.LPos<<endl; } } continue; } //txtFileNonPC<<entry<<endl; //filter out events with a low traceZ Float_t trace=999; Float_t traceZ=999; reco.CalcTrace(s,plInfo,&trace,&traceZ); if (traceZ<0.5){ cnt.badTraceZCounter++; //eventsTxtFileTrace<<event<<endl; MAXMSG("MeuAnalysis",Msg::kDebug,1000) <<"trace="<<trace<<", traceZ="<<traceZ <<", trk.EndTraceZ="<<trk.EndTraceZ<<endl; continue; } reco.CalcStripDists(s,plInfo); if (cuts.FilterBadDistEndStrip(s,plInfo)){//nothing for FD cnt.badDistEndCounter++; continue; } //MUST BE GOOD! cnt.goodWindowCounter++; cuts.ExtractTruthInfo(ev,s,plInfo); cuts.FillSTSumRecoDetails(s,plInfo); histos.FillMeuHistos(s); //eventsTxtFile<<entry<<endl; MAXMSG("MeuAnalysis",Msg::kInfo,10) <<"Good window: run=" <<ev.Run<<", snarl="<<ev.Snarl<<endl; if (s.MCLowEn>0.5*Munits::GeV) { //txtFileLowEn<<entry<<endl; MSG("MeuAnalysis",Msg::kInfo) <<"Low En="<<s.MCLowEn<<" high="<<s.MCHighEn <<", vtxPl="<<s.VtxPlane<<", endPl="<<s.EndPlane <<endl<<"entry="<<entry <<", snarl="<<s.Snarl<<", run="<<s.Run<<endl; } cuts.FillEnVsDist(s,plInfo); cuts.FillEnVsDist2(s,plInfo); //cuts.FillTimeHistos(ntp,evt,s); s.Temperature=this->GetTemperature(s); summaryOut.SummaryTreeFill(plInfo); }//end of for MSG("MeuAnalysis",Msg::kInfo)<<"Finished main loop"<<endl; MSG("MeuAnalysis",Msg::kInfo) <<endl<<"*** Run Summary ***"<<endl; histos.PrintMeuValues(); MSG("MeuAnalysis",Msg::kInfo) <<"Total Entries (snarls)="<<fEntries<<endl; cnt.PrintNtpSt(); //finish the summary tree and output it summaryOut.SummaryTreeFinish(); MSG("MeuAnalysis",Msg::kInfo) <<" ** Finished MakeSummaryTreeWithAtNu method **"<<endl; }

void MeuAnalysis::MakeSummaryTreeWithNtpSt (   )  [private]

Definition at line 5043 of file MeuAnalysis.cxx. References fRec, fRecBD, InitialiseLoopVariables(), MakeSummaryTreeWithNtpStOneSnarl(), MSG, SetLoopVariables(), and StoreOrFinishSummaryTree(). Referenced by MakeSummaryTree(). { MSG("MeuAnalysis",Msg::kInfo) <<" ** Running MakeSummaryTreeWithNtpSt method... **"<<endl; this->InitialiseLoopVariables(); //for(Int_t entry=55000;entry<fEntries;entry++){ for(Int_t entry=0;entry<fEntries;entry++){ this->SetLoopVariables(entry); //run the algorithm to calc meu values for events in snarl this->MakeSummaryTreeWithNtpStOneSnarl(fRec,fRecBD); }//end of for MSG("MeuAnalysis",Msg::kInfo)<<"Finished main loop"<<endl; //finish the job this->StoreOrFinishSummaryTree(NULL,NULL,1); MSG("MeuAnalysis",Msg::kInfo) <<" ** Finished MakeSummaryTreeWithNtpSt method **"<<endl; }

void MeuAnalysis::MakeSummaryTreeWithNtpStOneSnarl (  const NtpStRecord *  pntp, const NtpBDLiteRecord *  pntpBD )  const

code to produce the MeuSummary ntuple this should be called once per snarl Definition at line 4661 of file MeuAnalysis.cxx. References MeuCounter::bad2ndContCounter, MeuCounter::badDistEndCounter, MeuCounter::badEndGapsCounter, MeuCounter::badFidCounter, MeuCounter::badRecoCounter, MeuCounter::badShwDistCounter, MeuCounter::badTraceZCounter, MeuCounter::badTrackEndCounter, MeuCounter::badViewsCounter, MeuCounter::badWindowCounter, NtpSRPlane::beg, MeuCounter::bothSMHitCounter, MeuReco::CalcFCPC(), MeuReco::CalcStripDists(), MeuReco::CalcTrace(), MeuReco::CalcVtx(), MeuReco::CalcVtxSpecialND(), MeuReco::CalcWindow(), MeuCuts::CheckTrkViews(), MeuReco::CheckWinContainment(), MeuCounter::closeTimesCounter, MeuCounter::coilHitCounter, MeuCounter::correctTrigSrcCounter, MeuSummary::Detector, MeuSummary::DistToEdgeFid, NtpSRPlane::end, MeuSummary::EndDistToEdge, MeuSummary::EndPlane, MeuCounter::endPlaneSpectCounter, MeuCounter::entriesAll, MeuCounter::entry, MeuSummary::Event, NtpStRecord::evt, MeuCounter::evtCounter, MeuCuts::ExtractPlInfo(), MeuCuts::ExtractTruthInfo(), MeuSummary::FC, MeuCounter::FCCounter, MeuCuts::FillBeamMonDetails(), MeuCuts::FillEnVsDist(), MeuCuts::FillEnVsDist2(), MeuHistos::FillMeuHistos(), MeuCuts::FillSTSumDetails(), MeuCuts::FillSTSumRecoDetails(), MeuCuts::FillTimeHistos(), MeuCuts::FilterBadDistEndStrip(), MeuCuts::FilterBadEndGaps(), MeuCuts::FilterBadEvtPerSlc(), MeuCuts::FilterBadShwDist(), MeuCuts::FilterBadTrackEnd(), MeuCuts::FilterBadTrkTimes(), MeuCuts::GetBDSelectSpillInfo(), RecRecordImp< T >::GetHeader(), MeuSummaryWriter::GetMeuSummaryToFill(), RecDataHeader::GetRun(), RecPhysicsHeader::GetSnarl(), RecDataHeader::GetSubRun(), MeuCounter::goodContainmentCounter, MeuCounter::goodSpillCounter, MeuCounter::goodWindowCounter, MeuCuts::IsAwayFromCoil(), MeuCuts::IsCorrectTrigSrc(), MeuCounter::longTrack2Counter, MeuCounter::longTrackCounter, MeuHitInfo::LPos, MAXMSG, MeuSummary::MCHighEn, MeuSummary::MCLowEn, MSG, MeuCounter::notEmptyEventCounter, MeuCounter::notLICounter, NtpSRTrack::nstrip, NtpSREvent::ntrack, MeuCounter::oneEvtPerSlcCounter, MeuCounter::oneTrackCounter, MeuSummary::PC, MeuCounter::PCCounter, MeuHitInfo::Plane, NtpSRTrack::plane, print(), NtpStRecord::Print(), MeuCuts::PrintNtpSt(), MeuReco::Reconstruct(), MeuSummary::Run, run(), s(), MeuCuts::SetSpecificCuts(), NtpSREvent::slc, MeuSummary::SMBoth, MeuSummary::Snarl, MeuSummary::SubRun, MeuSummaryWriter::SummaryTreeFill(), MeuSummaryWriter::SummaryTreeSetup(), MeuSummary::Temperature, MeuCounter::TGCounter, MeuHitInfo::TPos, NtpSREvent::trk, NtpStRecord::trk, MeuSummary::VtxDistToEdge, MeuSummary::VtxPlane, MeuSummary::WinSigMap, MeuHitInfo::X, MeuHitInfo::Y, MeuHitInfo::Z, and MeuCounter::zeroEvtsCounter. Referenced by MeuCalModule::Ana(), and MakeSummaryTreeWithNtpSt(). { //get a reference (coding preference) //don't get a reference to NtpBDLiteRecord because it doesn't exist //for MC const NtpStRecord& ntp=*pntp; //get candheader const RecCandHeader& rec=ntp.GetHeader(); Int_t run=rec.GetRun(); Int_t subrun=rec.GetSubRun(); Int_t snarl=rec.GetSnarl(); static Int_t lastRun=-1; static Int_t lastSubrun=-1; static Int_t runningTotalOfRuns=0; static Int_t runningTotalOfSubruns=0; if (run!=lastRun) { runningTotalOfRuns++; MSG("MeuAnalysis",Msg::kInfo) <<"Found new run="<<run <<", running total of runs used="<<runningTotalOfRuns<<endl; } lastRun=run; if (subrun!=lastSubrun) { runningTotalOfSubruns++; MSG("MeuAnalysis",Msg::kInfo) <<"Found new subrun="<<subrun <<", running total of subruns used="<<runningTotalOfSubruns<<endl; } lastSubrun=subrun; MAXMSG("MeuAnalysis",Msg::kInfo,1) <<"First entry: run="<<run<<", snarl="<<snarl<<endl; static const Int_t planeCut=8; //variable to turn on all the useful messages if required static const Msg::LogLevel_t logLevel=Msg::kDebug; //create static variables to prevent multiple constructions static Bool_t firstTime=true; //reconstruction object static const MeuReco reco; static const MeuCuts cuts; static const MeuHistos histos; static MeuCounter cnt; //count all the calls of this loop cnt.entriesAll++; static MeuSummaryWriter* sOut=0; static ofstream* peventsTxtFile=0; static ofstream* ptxtFileLowEn=0; static ofstream* ptxtFileNonPC=0; if (firstTime) { //set the control variable so code not run again firstTime=false; //create the writer sOut=new MeuSummaryWriter(); sOut->SummaryTreeSetup(run,subrun); //this stores the objects' pointers for access at job end this->StoreOrFinishSummaryTree(sOut,&cnt,false); //printout the first element from the tree ntp.Print(cout); string temps="myeventsLowEn"; if (run>10000) temps+="MC"; peventsTxtFile=this->OpenTxtFile(100,"myevents"); ptxtFileLowEn=this->OpenTxtFile(100,temps.c_str()); ptxtFileNonPC=this->OpenTxtFile(100,"myeventsNonPC"); } //get references (coding preference) MeuSummaryWriter& summaryOut=*sOut; ofstream& eventsTxtFile=*peventsTxtFile; ofstream& txtFileLowEn=*ptxtFileLowEn; ofstream& txtFileNonPC=*ptxtFileNonPC; MSG("MeuAnalysis",logLevel) <<"Entry="<<cnt.entry<<endl; TClonesArray& evtTca=(*ntp.evt); const Int_t numEvts=evtTca.GetEntriesFast(); //count the number of snarls with zero events if (numEvts<=0) cnt.zeroEvtsCounter++; //loop over the events for (Int_t ievt=0;ievt<numEvts;++ievt) { const NtpSREvent* pevt= dynamic_cast<NtpSREvent*>(evtTca[ievt]); const NtpSREvent& evt=(*pevt); cnt.evtCounter++; //ensure there is only 1 track in the event if (evt.ntrack!=1) continue; cnt.oneTrackCounter++; //filter out events occuring in a slice with more than 1 event if (cuts.FilterBadEvtPerSlc(ntp,evt.slc,cnt.entry)) continue; cnt.oneEvtPerSlcCounter++; TClonesArray& trkTca=(*ntp.trk); const NtpSRTrack* ptrk= dynamic_cast<NtpSRTrack*>(trkTca[evt.trk[0]]); const NtpSRTrack& trk=(*ptrk); //cut on number of planes if (sqrt(pow(1.*trk.plane.beg- trk.plane.end,2))<planeCut) continue; cnt.longTrackCounter++; //MeuSummary s; MeuSummary& s=summaryOut.GetMeuSummaryToFill(); s.Event=cnt.entry; cuts.FillSTSumDetails(ntp,s,ievt,evt.slc); cuts.SetSpecificCuts(s); //map to hold the info for each plane in the entry map<Int_t,MeuHitInfo> plInfo; if (!cuts.IsCorrectTrigSrc(s)){ continue; } cnt.correctTrigSrcCounter++; //cut out the LI events if (this->FilterLI(ntp,s)){ //eventsTxtFileLI<<event<<endl; continue; } cnt.notLICounter++; if (pntpBD) cuts.GetBDSelectSpillInfo(*pntpBD,s,cnt.entry, cnt.goodSpillCounter); //print out info for specific snarl(s) if (s.Snarl==-1){ MSG("MeuAnalysis",Msg::kInfo) <<"Entry="<<cnt.entry<<", run="<<run<<", snarl="<<snarl <<", trk.nstrips="<<trk.nstrip<<endl; cuts.PrintNtpSt(ntp); } if (cuts.FilterBadTrkTimes(ntp,s,evt)) continue; cnt.closeTimesCounter++; cuts.ExtractPlInfo(ntp,s,plInfo,evt); if (plInfo.size()==0) { MAXMSG("MeuAnalysis",Msg::kWarning,20) <<"No data extracted for entry="<<cnt.entry<<endl; continue; } cnt.notEmptyEventCounter++; reco.CalcVtx(s,plInfo); //check if track stops in spectrometer if (s.Detector==Detector::kNear && s.EndPlane>120) continue; cnt.endPlaneSpectCounter++; MSG("MeuAnalysis",logLevel)<<"CalcFCPC..."<<endl; reco.CalcFCPC(s,plInfo); //recalculate the vtx of the tracks //coming in from the spectrometer MSG("MeuAnalysis",logLevel)<<"CalcVtxSpecialND..."<<endl; reco.CalcVtxSpecialND(s,plInfo); //cut on number of planes with potential new vertex if (sqrt(pow(1.*s.VtxPlane-s.EndPlane,2))<planeCut) { MAXMSG("MeuAnalysis",Msg::kDebug,100) <<"Cutting on track length, run="<<s.Run <<", subrun="<<s.SubRun<<", snarl="<<s.Snarl<<endl <<"s.VtxPlane="<<s.VtxPlane <<", s.EndPlane="<<s.EndPlane <<", trk.beg="<<trk.plane.beg <<", trk.end="<<trk.plane.end<<endl; //cuts.PrintNtpSt(ntp); continue; } cnt.longTrack2Counter++; Bool_t TG=!s.FC && !s.PC; if (s.PC) cnt.PCCounter++; else if (s.FC) cnt.FCCounter++; else if (TG) cnt.TGCounter++; if (s.FC && s.PC) MSG("MeuAnalysis",Msg::kError) <<"Both FC and PC!"<<endl; Bool_t awayFromCoil=cuts.IsAwayFromCoil(s,plInfo); if (s.PC && !awayFromCoil) cnt.coilHitCounter++; if (s.PC && s.SMBoth && awayFromCoil) cnt.bothSMHitCounter++; Bool_t goodContainment=s.PC && !s.SMBoth && awayFromCoil; if (!goodContainment) continue; cnt.goodContainmentCounter++; MSG("MeuAnalysis",logLevel)<<"Found Good PC event..."<<endl; Int_t diffAtEnd=-1; Int_t diffAtVtx=-1; if (!cuts.CheckTrkViews(s,&diffAtVtx,&diffAtEnd)) { cnt.badViewsCounter++; continue; } //check that the hits beyond the end of the track are ok if (cuts.FilterBadTrackEnd(s,plInfo)){ cnt.badTrackEndCounter++; continue; } Bool_t goodReco=reco.Reconstruct(s,plInfo); if (!goodReco){ cnt.badRecoCounter++; MSG("MeuAnalysis",logLevel)<<"Bad reco..."<<endl; continue; } //check that there are no big gaps at the end of the track if (cuts.FilterBadEndGaps(s,plInfo)){//after reco only cnt.badEndGapsCounter++; continue; } //recalculate the containment (using MeuReco reconstruction) s.DistToEdgeFid=0.35;//allow a slight shift in position reco.CalcFCPC(s,plInfo); //re-check containment (using MeuReco reconstruction) if (!s.PC){ cnt.bad2ndContCounter++; static Bool_t print=false; if (print){ MAXMSG("MeuAnalysis",Msg::kInfo,100) <<endl<<endl<<"Recalculated containment and not PC!" <<" PC="<<s.PC<<", FC="<<s.FC <<" entry="<<cnt.entry <<", snarl="<<s.Snarl<<", run="<<s.Run<<endl <<"Low En="<<s.MCLowEn<<" high="<<s.MCHighEn<<endl <<endl; MAXMSG("MeuAnalysis",Msg::kInfo,100) <<"vtx Pl="<<s.VtxPlane<<", end Pl="<<s.EndPlane <<", dVtx="<<s.VtxDistToEdge <<", dEnd="<<s.EndDistToEdge<<endl; for (map<Int_t,MeuHitInfo>::const_iterator it= plInfo.begin(); it!=plInfo.end();++it){ const MeuHitInfo& c=it->second; MAXMSG("MeuAnalysis",Msg::kInfo,500) <<"pl="<<c.Plane<<", z="<<c.Z <<", x="<<c.X<<", y="<<c.Y <<" t="<<c.TPos<<", l="<<c.LPos<<endl; } } txtFileNonPC<<cnt.entry<<endl; continue; } MSG("MeuAnalysis",logLevel)<<"Recalibrating..."<<endl; this->Recalibrate(ntp,s,plInfo,evt); MSG("MeuAnalysis",logLevel)<<"CalcWindow..."<<endl; reco.CalcWindow(s,plInfo); MSG("MeuAnalysis",logLevel)<<"CheckContainment..."<<endl; //check that a track window was found if (s.WinSigMap<=0) { cnt.badWindowCounter++; continue; } //check that the track window is contained within the fid vol if (!reco.CheckWinContainment(s,plInfo)){//nothing for FD cnt.badFidCounter++; continue; } //filter out events with a low traceZ Float_t trace=999; Float_t traceZ=999; reco.CalcTrace(s,plInfo,&trace,&traceZ);//only FD now if (traceZ<0.5){ cnt.badTraceZCounter++; //eventsTxtFileTrace<<event<<endl; MAXMSG("MeuAnalysis",Msg::kDebug,1000) <<"trace="<<trace<<", traceZ="<<traceZ<<endl; continue; } //check that all the hits in the track window are >20 cm from //the end of the strips reco.CalcStripDists(s,plInfo); if (cuts.FilterBadDistEndStrip(s,plInfo)){//requires window cnt.badDistEndCounter++; continue; } //if beam event check that the track window is away from //the possible hadronic shower if (cuts.FilterBadShwDist(ntp,s,evt)){//only for ND spill data cnt.badShwDistCounter++; continue; } //MUST BE GOOD! cnt.goodWindowCounter++; cuts.ExtractTruthInfo(ntp,s,plInfo); cuts.FillSTSumRecoDetails(s,plInfo); histos.FillMeuHistos(s); if (pntpBD) cuts.FillBeamMonDetails(*pntpBD,s); eventsTxtFile<<cnt.entry<<endl; MAXMSG("MeuAnalysis",Msg::kInfo,10) <<"Good track window: run="<<run<<", snarl="<<snarl <<", evt="<<ievt+1<<"/"<<numEvts<<endl; if (s.MCLowEn>0.5*Munits::GeV) { txtFileLowEn<<cnt.entry<<endl; MSG("MeuAnalysis",Msg::kInfo) <<"Low En="<<s.MCLowEn<<" high="<<s.MCHighEn <<", vtxPl="<<s.VtxPlane<<", endPl="<<s.EndPlane <<endl<<"entry="<<cnt.entry <<", snarl="<<s.Snarl<<", run="<<s.Run<<endl; } cuts.FillEnVsDist(s,plInfo); cuts.FillEnVsDist2(s,plInfo); cuts.FillTimeHistos(ntp,evt,s);//a loop over lots of strips... s.Temperature=this->GetTemperature(s); summaryOut.SummaryTreeFill(plInfo); }//end of loop over events //increment ntuple entry (snarl) counter cnt.entry++; }

void MeuAnalysis::N_1Plots (   )

Definition at line 1114 of file MeuAnalysis.cxx. References abs(), MeuReco::CalcFCPC(), MeuReco::CalcSmallestDeepDistToEdge(), MeuReco::CalcStripDists(), MeuReco::CalcVtx(), MeuReco::CalcVtxSpecialND(), MeuReco::CalcWindow(), MeuCuts::CheckTrkViews(), MeuReco::CheckWinContainment(), MeuSummary::DistToEdgeFid, NtpSRTrack::end, MeuSummary::EndPlane, MeuSummary::Event, NtpStRecord::evt, MeuCuts::ExtractPlInfo(), MeuCuts::FillSTSumDetails(), MeuCuts::FilterBadDistEndStrip(), MeuCuts::FilterBadEndGaps(), MeuCuts::FilterBadTrackEnd(), MeuCuts::FilterLowMatTraversed(), fOutFile, InitialiseLoopVariables(), MAXMSG, MSG, OpenFile(), MeuSummary::PC, NtpSRVertex::plane, Recalibrate(), MeuReco::Reconstruct(), MeuSummary::Reset(), MeuSummary::RFid, s(), SetLoopVariables(), MeuSummary::TotalMatTraversed, NtpStRecord::trk, NtpSRTrack::vtx, MeuSummary::VtxPlane, MeuSummary::WinAvPLCor, MeuSummary::WinSigMap, MeuSummary::WinStopSidePl, and MeuSummary::WinVtxSidePl. { MSG("MeuAnalysis",Msg::kInfo) <<" ** Running N_1Plots method... **"<<endl; //open the output file for the histograms fOutFile=this->OpenFile(100,"N_1Plots"); TH1F* hTrkPl=new TH1F("hTrkPl","hTrkPl",201,-1,200); hTrkPl->SetTitle("hTrkPl"); hTrkPl->GetXaxis()->SetTitle("hTrkPl"); hTrkPl->GetXaxis()->CenterTitle(); hTrkPl->SetFillColor(0); hTrkPl->SetLineWidth(3); //hTrkPl->SetBit(TH1::kCanRebin); TH1F* hTrkPlN_1=new TH1F("hTrkPlN_1","hTrkPlN_1",201,-1,200); hTrkPlN_1->SetTitle("hTrkPlN_1"); hTrkPlN_1->GetXaxis()->SetTitle("hTrkPlN_1"); hTrkPlN_1->GetXaxis()->CenterTitle(); hTrkPlN_1->SetLineColor(2);//red for N-1 cuts hTrkPlN_1->SetFillColor(0); hTrkPlN_1->SetLineWidth(3); //hTrkPlN_1->SetBit(TH1::kCanRebin); TH1F* hTrkPlN=new TH1F("hTrkPlN","hTrkPlN",201,-1,200); hTrkPlN->SetTitle("hTrkPlN"); hTrkPlN->GetXaxis()->SetTitle("hTrkPlN"); hTrkPlN->GetXaxis()->CenterTitle(); hTrkPlN->SetLineColor(4);//blue for N cuts hTrkPlN->SetFillColor(0); hTrkPlN->SetLineWidth(3); //hTrkPlN->SetBit(TH1::kCanRebin); TH1F* hDistEdge=new TH1F("hDistEdge","hDistEdge",201,-3,3); hDistEdge->SetTitle("hDistEdge"); hDistEdge->GetXaxis()->SetTitle("hDistEdge"); hDistEdge->GetXaxis()->CenterTitle(); hDistEdge->SetFillColor(0); hDistEdge->SetLineWidth(3); //hDistEdge->SetBit(TH1::kCanRebin); TH1F* hDistEdgeN_1=new TH1F("hDistEdgeN_1","hDistEdgeN_1",201,-3,3); hDistEdgeN_1->SetTitle("hDistEdgeN_1"); hDistEdgeN_1->GetXaxis()->SetTitle("hDistEdgeN_1"); hDistEdgeN_1->GetXaxis()->CenterTitle(); hDistEdgeN_1->SetLineColor(2);//red for N-1 cuts hDistEdgeN_1->SetFillColor(0); hDistEdgeN_1->SetLineWidth(3); //hDistEdgeN_1->SetBit(TH1::kCanRebin); TH1F* hDistEdgeN=new TH1F("hDistEdgeN","hDistEdgeN",201,-3,3); hDistEdgeN->SetTitle("hDistEdgeN"); hDistEdgeN->GetXaxis()->SetTitle("hDistEdgeN"); hDistEdgeN->GetXaxis()->CenterTitle(); hDistEdgeN->SetLineColor(4);//blue for N cuts hDistEdgeN->SetFillColor(0); hDistEdgeN->SetLineWidth(3); //hDistEdgeN->SetBit(TH1::kCanRebin); TH1F* hViewsVtx=new TH1F("hViewsVtx","hViewsVtx",101,-1,100); hViewsVtx->SetTitle("hViewsVtx"); hViewsVtx->GetXaxis()->SetTitle("hViewsVtx"); hViewsVtx->GetXaxis()->CenterTitle(); hViewsVtx->SetFillColor(0); hViewsVtx->SetLineWidth(3); //hViewsVtx->SetBit(TH1::kCanRebin); TH1F* hViewsVtxN_1=new TH1F("hViewsVtxN_1","hViewsVtxN_1",101,-1,100); hViewsVtxN_1->SetTitle("hViewsVtxN_1"); hViewsVtxN_1->GetXaxis()->SetTitle("hViewsVtxN_1"); hViewsVtxN_1->GetXaxis()->CenterTitle(); hViewsVtxN_1->SetLineColor(2);//red for N-1 cuts hViewsVtxN_1->SetFillColor(0); hViewsVtxN_1->SetLineWidth(3); //hViewsVtxN_1->SetBit(TH1::kCanRebin); TH1F* hViewsVtxN=new TH1F("hViewsVtxN","hViewsVtxN",101,-1,100); hViewsVtxN->SetTitle("hViewsVtxN"); hViewsVtxN->GetXaxis()->SetTitle("hViewsVtxN"); hViewsVtxN->GetXaxis()->CenterTitle(); hViewsVtxN->SetLineColor(4);//blue for N cuts hViewsVtxN->SetFillColor(0); hViewsVtxN->SetLineWidth(3); //hViewsVtxN->SetBit(TH1::kCanRebin); TH1F* hViewsEnd=new TH1F("hViewsEnd","hViewsEnd",101,-1,100); hViewsEnd->SetTitle("hViewsEnd"); hViewsEnd->GetXaxis()->SetTitle("hViewsEnd"); hViewsEnd->GetXaxis()->CenterTitle(); hViewsEnd->SetFillColor(0); hViewsEnd->SetLineWidth(3); //hViewsEnd->SetBit(TH1::kCanRebin); TH1F* hViewsEndN_1=new TH1F("hViewsEndN_1","hViewsEndN_1",101,-1,100); hViewsEndN_1->SetTitle("hViewsEndN_1"); hViewsEndN_1->GetXaxis()->SetTitle("hViewsEndN_1"); hViewsEndN_1->GetXaxis()->CenterTitle(); hViewsEndN_1->SetLineColor(2);//red for N-1 cuts hViewsEndN_1->SetFillColor(0); hViewsEndN_1->SetLineWidth(3); //hViewsEndN_1->SetBit(TH1::kCanRebin); TH1F* hViewsEndN=new TH1F("hViewsEndN","hViewsEndN",101,-1,100); hViewsEndN->SetTitle("hViewsEndN"); hViewsEndN->GetXaxis()->SetTitle("hViewsEndN"); hViewsEndN->GetXaxis()->CenterTitle(); hViewsEndN->SetLineColor(4);//blue for N cuts hViewsEndN->SetFillColor(0); hViewsEndN->SetLineWidth(3); //hViewsEndN->SetBit(TH1::kCanRebin); TH1F* hDistEndSt=new TH1F("hDistEndSt","hDistEndSt",200,-0.1,2); hDistEndSt->SetTitle("hDistEndSt"); hDistEndSt->GetXaxis()->SetTitle("hDistEndSt"); hDistEndSt->GetXaxis()->CenterTitle(); hDistEndSt->SetFillColor(0); hDistEndSt->SetLineWidth(3); //hDistEndSt->SetBit(TH1::kCanRebin); TH1F* hDistEndStN_1=new TH1F("hDistEndStN_1","hDistEndStN_1", 200,-0.1,2); hDistEndStN_1->SetTitle("hDistEndStN_1"); hDistEndStN_1->GetXaxis()->SetTitle("hDistEndStN_1"); hDistEndStN_1->GetXaxis()->CenterTitle(); hDistEndStN_1->SetLineColor(2);//red for N-1 cuts hDistEndStN_1->SetFillColor(0); hDistEndStN_1->SetLineWidth(3); //hDistEndStN_1->SetBit(TH1::kCanRebin); TH1F* hDistEndStN=new TH1F("hDistEndStN","hDistEndStN",200,-0.1,2); hDistEndStN->SetTitle("hDistEndStN"); hDistEndStN->GetXaxis()->SetTitle("hDistEndStN"); hDistEndStN->GetXaxis()->CenterTitle(); hDistEndStN->SetLineColor(4);//blue for N cuts hDistEndStN->SetFillColor(0); hDistEndStN->SetLineWidth(3); //hDistEndStN->SetBit(TH1::kCanRebin); TH1F* hSigBeyondEnd=new TH1F("hSigBeyondEnd","hSigBeyondEnd",1000,-1,50000); hSigBeyondEnd->SetTitle("hSigBeyondEnd"); hSigBeyondEnd->GetXaxis()->SetTitle("hSigBeyondEnd"); hSigBeyondEnd->GetXaxis()->CenterTitle(); hSigBeyondEnd->SetFillColor(0); hSigBeyondEnd->SetLineWidth(3); //hSigBeyondEnd->SetBit(TH1::kCanRebin); TH1F* hSigBeyondEndN_1=new TH1F("hSigBeyondEndN_1","hSigBeyondEndN_1", 1000,-1,50000); hSigBeyondEndN_1->SetTitle("hSigBeyondEndN_1"); hSigBeyondEndN_1->GetXaxis()->SetTitle("hSigBeyondEndN_1"); hSigBeyondEndN_1->GetXaxis()->CenterTitle(); hSigBeyondEndN_1->SetLineColor(2);//red for N-1 cuts hSigBeyondEndN_1->SetFillColor(0); hSigBeyondEndN_1->SetLineWidth(3); //hSigBeyondEndN_1->SetBit(TH1::kCanRebin); TH1F* hSigBeyondEndN=new TH1F("hSigBeyondEndN","hSigBeyondEndN", 1000,-1,50000); hSigBeyondEndN->SetTitle("hSigBeyondEndN"); hSigBeyondEndN->GetXaxis()->SetTitle("hSigBeyondEndN"); hSigBeyondEndN->GetXaxis()->CenterTitle(); hSigBeyondEndN->SetLineColor(4);//blue for N cuts hSigBeyondEndN->SetFillColor(0); hSigBeyondEndN->SetLineWidth(3); //hSigBeyondEndN->SetBit(TH1::kCanRebin); TH1F* hSigTrkEndGaps=new TH1F("hSigTrkEndGaps","hSigTrkEndGaps",1000,-1,50000); hSigTrkEndGaps->SetTitle("hSigTrkEndGaps"); hSigTrkEndGaps->GetXaxis()->SetTitle("hSigTrkEndGaps"); hSigTrkEndGaps->GetXaxis()->CenterTitle(); hSigTrkEndGaps->SetFillColor(0); hSigTrkEndGaps->SetLineWidth(3); //hSigTrkEndGaps->SetBit(TH1::kCanRebin); TH1F* hSigTrkEndGapsN_1=new TH1F("hSigTrkEndGapsN_1","hSigTrkEndGapsN_1", 1000,-1,50000); hSigTrkEndGapsN_1->SetTitle("hSigTrkEndGapsN_1"); hSigTrkEndGapsN_1->GetXaxis()->SetTitle("hSigTrkEndGapsN_1"); hSigTrkEndGapsN_1->GetXaxis()->CenterTitle(); hSigTrkEndGapsN_1->SetLineColor(2);//red for N-1 cuts hSigTrkEndGapsN_1->SetFillColor(0); hSigTrkEndGapsN_1->SetLineWidth(3); //hSigTrkEndGapsN_1->SetBit(TH1::kCanRebin); TH1F* hSigTrkEndGapsN=new TH1F("hSigTrkEndGapsN","hSigTrkEndGapsN", 1000,-1,50000); hSigTrkEndGapsN->SetTitle("hSigTrkEndGapsN"); hSigTrkEndGapsN->GetXaxis()->SetTitle("hSigTrkEndGapsN"); hSigTrkEndGapsN->GetXaxis()->CenterTitle(); hSigTrkEndGapsN->SetLineColor(4);//blue for N cuts hSigTrkEndGapsN->SetFillColor(0); hSigTrkEndGapsN->SetLineWidth(3); //hSigTrkEndGapsN->SetBit(TH1::kCanRebin); TH1F* hMatTrav=new TH1F("hMatTrav","hMatTrav",400,0,10); hMatTrav->SetTitle("hMatTrav"); hMatTrav->GetXaxis()->SetTitle("hMatTrav"); hMatTrav->GetXaxis()->CenterTitle(); hMatTrav->SetFillColor(0); hMatTrav->SetLineWidth(3); //hMatTrav->SetBit(TH1::kCanRebin); TH1F* hMatTravN_1=new TH1F("hMatTravN_1","hMatTravN_1",400,0,10); hMatTravN_1->SetTitle("hMatTravN_1"); hMatTravN_1->GetXaxis()->SetTitle("hMatTravN_1"); hMatTravN_1->GetXaxis()->CenterTitle(); hMatTravN_1->SetLineColor(2);//red for N-1 cuts hMatTravN_1->SetFillColor(0); hMatTravN_1->SetLineWidth(3); //hMatTravN_1->SetBit(TH1::kCanRebin); TH1F* hMatTravN=new TH1F("hMatTravN","hMatTravN",400,0,10); hMatTravN->SetTitle("hMatTravN"); hMatTravN->GetXaxis()->SetTitle("hMatTravN"); hMatTravN->GetXaxis()->CenterTitle(); hMatTravN->SetLineColor(4);//blue for N cuts hMatTravN->SetFillColor(0); hMatTravN->SetLineWidth(3); //hMatTravN->SetBit(TH1::kCanRebin); //reconstruction object MeuReco reco; MeuCuts cuts; MeuSummary s; this->InitialiseLoopVariables(); //for(Int_t entry=0;entry<10000;entry++){ for(Int_t entry=0;entry<fEntries;entry++){ this->SetLoopVariables(entry); NtpStRecord& ntp=(*fRec);//Make a reference instead of a pointer //ensure there is only 1 track in the entry TClonesArray& trkTca=(*ntp.trk); Int_t numTrks=trkTca.GetEntries(); if (numTrks!=1) continue; const NtpSRTrack* ptrk=dynamic_cast<NtpSRTrack*>(trkTca[0]); const NtpSRTrack& trk=(*ptrk); //cut on number of planes Int_t numPlanes=abs(trk.vtx.plane-trk.end.plane); if (numPlanes<8) continue; s.Reset(); s.RFid=1.0; s.DistToEdgeFid=0.4; s.Event=entry; cuts.FillSTSumDetails(ntp,s,-1,-1); map<Int_t,MeuHitInfo> plInfo; //check if the track is in the spectrometer if (trk.vtx.plane>120 && trk.end.plane>120) continue; TClonesArray& evtTca=(*ntp.evt); const NtpSREvent* pevt=dynamic_cast<NtpSREvent*>(evtTca[0]); cuts.ExtractPlInfo(ntp,s,plInfo,*pevt); reco.CalcVtx(s,plInfo); //check if track stops in spectrometer if (s.EndPlane>120) continue; //calculate the containment (before doing the vtxSpecial) reco.CalcFCPC(s,plInfo); //recalculate the vtx of the tracks coming in from the spectrometer reco.CalcVtxSpecialND(s,plInfo); //cut on number of planes with potential new vertex numPlanes=abs(s.VtxPlane-s.EndPlane); if (numPlanes<8) continue; reco.CalcStripDists(s,plInfo); //calc distance to detector edge Float_t dist=reco.CalcSmallestDeepDistToEdge(plInfo[s.EndPlane]); Int_t diffAtEnd=-1; Int_t diffAtVtx=-1; cuts.CheckTrkViews(s,&diffAtVtx,&diffAtEnd); Float_t minDist=-1; Float_t sigBeyondEnd=-1; Float_t adcInGap=-1; Bool_t goodContainment=false; Bool_t goodReco=false; if (s.PC) { goodReco=reco.Reconstruct(s,plInfo); if (goodReco) { TClonesArray& evtTca=(*ntp.evt); const NtpSREvent* pevt=dynamic_cast<NtpSREvent*>(evtTca[0]); this->Recalibrate(ntp,s,plInfo,*pevt); reco.CalcWindow(s,plInfo); //now check that window is deeply contained if (s.WinSigMap>0) { goodContainment=reco.CheckWinContainment(s,plInfo); } //now do some of the calculations that require good PC if (goodContainment){ cuts.FilterBadDistEndStrip(s,plInfo,&minDist); cuts.FilterBadTrackEnd(s,plInfo,&sigBeyondEnd); cuts.FilterBadEndGaps(s,plInfo,&adcInGap); } } } //if (s.WinSigMap<=0) { //matTrav //start of "no cuts" hTrkPl->Fill(numPlanes); hDistEdge->Fill(dist); hViewsVtx->Fill(diffAtVtx); hViewsEnd->Fill(diffAtEnd); if (goodContainment){ //only have a value in this case hDistEndSt->Fill(minDist); hSigBeyondEnd->Fill(sigBeyondEnd); hSigTrkEndGaps->Fill(adcInGap); } //only have a value in this case (requires calcWindow) if (goodReco) hMatTrav->Fill(s.TotalMatTraversed); //cuts Bool_t goodDist=dist>0.4; Bool_t goodTrkPl=numPlanes>10; Bool_t goodViewsVtx=diffAtVtx<=9;//for ND full inst. planes Bool_t goodViewsEnd=diffAtEnd<=5; Bool_t goodDistEndSt=minDist>0.2; Bool_t goodSigBeyondEnd=sigBeyondEnd<500 && sigBeyondEnd>=0; Bool_t goodSigTrkEndGaps=adcInGap<500 && adcInGap>=0; Bool_t goodMatTrav=!cuts.FilterLowMatTraversed(s); Bool_t goodN= goodContainment && goodDist && goodTrkPl && goodViewsVtx && goodViewsEnd && goodDistEndSt && goodSigBeyondEnd && goodSigTrkEndGaps && goodMatTrav && true;//for coding convenience //start of N-1 cuts plots if (goodContainment && goodDist && goodViewsVtx && goodViewsEnd && goodDistEndSt && goodSigBeyondEnd && goodSigTrkEndGaps && goodMatTrav && true) { hTrkPlN_1->Fill(numPlanes); if (numPlanes<10){ MAXMSG("MeuAnalysis",Msg::kInfo,10) <<"vtxPl="<<s.VtxPlane <<", endPl="<<s.EndPlane <<", stopPl="<<s.WinVtxSidePl <<", startPl="<<s.WinStopSidePl <<", meu="<<s.WinSigMap <<", PLC="<<s.WinAvPLCor <<endl; } } if (goodContainment && goodTrkPl && goodViewsVtx && goodViewsEnd && goodDistEndSt && goodSigBeyondEnd && goodSigTrkEndGaps && goodMatTrav && true) { hDistEdgeN_1->Fill(dist); } if (goodContainment && goodDist && goodTrkPl && goodViewsEnd && goodDistEndSt && goodSigBeyondEnd && goodSigTrkEndGaps && goodMatTrav && true) { hViewsVtxN_1->Fill(diffAtVtx); } if (goodContainment && goodDist && goodTrkPl && goodViewsVtx && goodDistEndSt && goodSigBeyondEnd && goodSigTrkEndGaps && goodMatTrav && true) { hViewsEndN_1->Fill(diffAtEnd); } if (goodContainment && goodDist && goodTrkPl && goodViewsVtx && goodViewsEnd && goodSigBeyondEnd && goodSigTrkEndGaps && goodMatTrav && true) { hDistEndStN_1->Fill(minDist); } if (goodContainment && goodDist && goodTrkPl && goodViewsVtx && goodViewsEnd && goodDistEndSt && goodSigTrkEndGaps && goodMatTrav && true) { hSigBeyondEndN_1->Fill(sigBeyondEnd); } if (goodContainment && goodDist && goodTrkPl && goodViewsVtx && goodViewsEnd && goodDistEndSt && goodSigBeyondEnd && goodMatTrav && true) { hSigTrkEndGapsN_1->Fill(adcInGap); } if (goodContainment && goodDist && goodTrkPl && goodViewsVtx && goodViewsEnd && goodDistEndSt && goodSigBeyondEnd && goodSigTrkEndGaps && true) { hMatTravN_1->Fill(s.TotalMatTraversed); } //start of N cuts plots if (goodN) { hTrkPlN->Fill(numPlanes); hDistEdgeN->Fill(dist); hViewsVtxN->Fill(diffAtVtx); hViewsEndN->Fill(diffAtEnd); hDistEndStN->Fill(minDist); hSigBeyondEndN->Fill(sigBeyondEnd); hSigTrkEndGapsN->Fill(adcInGap); hMatTravN->Fill(s.TotalMatTraversed); } }//end of for MSG("MeuAnalysis",Msg::kInfo)<<"Finished main loop"<<endl; MSG("MeuAnalysis",Msg::kInfo) <<" ** Finished N_1Plots method **"<<endl; }

Bool_t MeuAnalysis::ObjectExistsInFile (  TFile *  file, std::string  objectName )  const

Definition at line 338 of file MeuAnalysis.cxx. References MSG, and RecCandHeader::Print(). Referenced by MakeChain(). { if (!file) { MSG("MeuAnalysis",Msg::kWarning) <<"File does not exist so can't test if object exists!"<<endl; return false; } TList* listOfKeys=file->GetListOfKeys(); static TFile* lastFile=0; if (lastFile!=file){ MSG("MeuAnalysis",Msg::kInfo) <<"File contains these keys:"<<endl; listOfKeys->Print(); } lastFile=file; TObject* ob=listOfKeys->FindObject(objectName.c_str()); if (ob) { MSG("MeuAnalysis",Msg::kInfo) <<"Requested object with name="<<objectName <<" exists in file:"<<endl; ob->Print(); return true; } else { MSG("MeuAnalysis",Msg::kInfo) <<"Requested object with name="<<objectName <<" does NOT exist in file"<<endl; return false; } }

TFile * MeuAnalysis::OpenFile (  Int_t  runNumber, std::string  prefix )  const [private]

Definition at line 650 of file MeuAnalysis.cxx. References Form(), MSG, and Test(). Referenced by BasicPlots(), N_1Plots(), and SpillPlots(). { //create the tfile pointer to be returned TFile* outputFile=0; //get the environmental variable char *anaDir=getenv("MEUANA_DIR"); //use a string to hold env instead string sAnaDir=""; if (anaDir!=NULL) { sAnaDir=anaDir; } else { MSG("MeuAnalysis",Msg::kInfo) <<"Environmental variable $MEUANA_DIR not set." <<" Writing file(s) to current directory"<<endl; sAnaDir="."; } //convert variables to string string sRunNumber=Form("%d",runNumber); //string sDetector="C"; string sDetector=""; //string sPrefix="h";//default string sPrefix="";//default if (prefix!="") sPrefix+=prefix; string sBase=sAnaDir+"/"+sPrefix+sDetector+sRunNumber; string sFileName=sBase+".root"; //test if file already exists ifstream Test(sFileName.c_str()); //open the appropriate file if(!Test){ outputFile=new TFile(sFileName.c_str(),"RECREATE"); } else { //Need new filename Int_t fred=1; while(Test) { Test.close(); string sAppendage=Form("%d",fred); sFileName=sBase+"_"+sAppendage+".root"; Test.open(sFileName.c_str()); fred++; } outputFile=new TFile(sFileName.c_str(),"NEW"); outputFile->SetCompressionLevel(9); } string sTmp="No File!"; if (outputFile) sTmp=outputFile->GetName(); MSG("MeuAnalysis",Msg::kInfo) <<"Output file opened: "<<sTmp<<endl; return outputFile; }

ofstream * MeuAnalysis::OpenTxtFile (  Int_t  runNumber, std::string  prefix )  const [private]

Definition at line 712 of file MeuAnalysis.cxx. References Form(), and MSG. Referenced by BasicReco(), SnarlList(), and SpillPlots(). { //create the tfile pointer to be returned ofstream* outputFile=0; //convert variables to string string sRunNumber=Form("%d",runNumber); string sDetector="F"; string sPrefix="";//default string sAnaDir="."; if (prefix!="") sPrefix+=prefix; string sBase=sAnaDir+"/"+sPrefix+sDetector+sRunNumber; string sFileName=sBase+".txt"; outputFile=new ofstream(sFileName.c_str()); if (outputFile){ MSG("MeuAnalysis",Msg::kInfo) <<"Output txt file opened: "<<sFileName<<endl; } else{ MSG("MeuAnalysis",Msg::kInfo) <<"Txt file NOT opened: "<<sFileName<<endl; } return outputFile; }

Bool_t MeuAnalysis::PlaneIsWithinTrack (  const MeuSummary &  s, Int_t  plane, Int_t  vtxPl, Int_t  endPl )  const [private]

Definition at line 2046 of file MeuAnalysis.cxx. References MeuSummary::Detector, MSG, and s(). Referenced by Recalibrate(). { Bool_t withinTrk=false; if (s.Detector==Detector::kNear){ if (plane<0 || plane>120) { MSG("MeuAnalysis",Msg::kVerbose) <<"MeuAnalysis::PlaneIsWithinTrack." <<" Plane is out of range="<<plane <<", det="<<s.Detector<<endl; return withinTrk; } } else if (s.Detector==Detector::kFar){ if (plane<0 || plane>485) { MSG("MeuAnalysis",Msg::kVerbose) <<"MeuAnalysis::PlaneIsWithinTrack." <<" Plane is out of range="<<plane <<", det="<<s.Detector<<endl; return withinTrk; } } else { MSG("MeuAnalysis",Msg::kWarning) <<"MeuAnalysis::PlaneIsWithinTrack." <<" Detector not supported"<<plane <<endl; } //count the end planes as being "within" the track if ((plane>=vtxPl && plane<=endPl) || (plane<=vtxPl && plane>=endPl)) withinTrk=true; return withinTrk; }

void MeuAnalysis::Recalibrate (  const NtpStRecord &  ntp, MeuSummary &  s, std::map< Int_t, MeuHitInfo > &  plInfo, const NtpSREvent &  evt )  const [private]

always have to calibrate sigmap since the information is not available for strips that aren't in tracks or showers Definition at line 2126 of file MeuAnalysis.cxx. References MeuSummary::Detector, MeuSummary::EndPlane, Calibrator::GetAttenCorrectedTpos(), Calibrator::GetDriftCorrected(), Calibrator::GetLinearized(), Calibrator::GetPhotoElectrons(), Calibrator::GetStripToStripCorrected(), Calibrator::Instance(), MeuHitInfo::LPos, MAXMSG, MeuHitInfo::MCLPos, MeuHitInfo::MCSigMap, MeuSummary::MedianTime, NtpSRSlice::nstrip, NtpSRPulseHeight::pe, MeuHitInfo::Pe, MeuHitInfo::Pe1, MeuHitInfo::Pe2, NtpSRStrip::ph0, NtpSRStrip::ph1, NtpSRStrip::plane, MeuHitInfo::Plane, PlaneIsWithinTrack(), print(), NtpSRPulseHeight::raw, Calibrator::ReInitialise(), s(), NtpSRPulseHeight::sigcor, MeuHitInfo::SigCor, MeuHitInfo::SigCor1, MeuHitInfo::SigCor2, MeuHitInfo::SigCorTrk1, MeuHitInfo::SigCorTrk2, MeuHitInfo::SigDrf, MeuHitInfo::SigDrf1, MeuHitInfo::SigDrf2, NtpSRPulseHeight::siglin, MeuHitInfo::SigLin, MeuHitInfo::SigLin1, MeuHitInfo::SigLin2, MeuHitInfo::SigLinOnly, MeuHitInfo::SigLinOnly1, MeuHitInfo::SigLinOnly2, MeuHitInfo::SigMap, MeuHitInfo::SigMap1, MeuHitInfo::SigMap2, MeuSummary::SimFlag, NtpSREvent::slc, NtpStRecord::slc, NtpSRSlice::stp, NtpStRecord::stp, NtpSRStrip::strip, MeuHitInfo::Strip, NtpSRStrip::time0, NtpSRStrip::time1, MeuSummary::TimeSec, and MeuSummary::VtxPlane. Referenced by BasicReco(), N_1Plots(), and SpillPlots(). { //variable to turn on all the useful messages if required Msg::LogLevel_t logLevel=Msg::kDebug; //Msg::LogLevel_t logLevelSts=Msg::kDebug; //Msg::LogLevel_t logLevelLin=Msg::kInfo; //make plots of the calibration constants for the recalibration //this is also done elsewhere for ntuple values of the cal constants //if the recalibration is not done then they are the same as sntp static TH1F* hAdcPerSigLinReCal1=0; static TH1F* hSigLinPerSigCorReCal1=0; static TH1F* hSigCorPerSigMapReCal1=0; static TH1F* hSigMapPerMipReCal1=0; static TH1F* hMipPerGeVReCal1=0; static TH1F* hAdcPerPeReCal1=0; static TH1F* hAdcPerSigDrfReCal1=0; static TH1F* hAdcPerSigLinOnlyReCal1=0; static TH1F* hAdcPerSigLinReCal2=0; static TH1F* hSigLinPerSigCorReCal2=0; static TH1F* hSigCorPerSigMapReCal2=0; static TH1F* hSigMapPerMipReCal2=0; static TH1F* hMipPerGeVReCal2=0; static TH1F* hAdcPerPeReCal2=0; static TH1F* hAdcPerSigDrfReCal2=0; static TH1F* hAdcPerSigLinOnlyReCal2=0; if (hAdcPerSigLinReCal1==0) { //stripend 1 hAdcPerSigLinReCal1=new TH1F("hAdcPerSigLinReCal1", "hAdcPerSigLinReCal1", 1100,-1,10); hAdcPerSigLinReCal1->SetFillColor(0); hAdcPerSigLinReCal1->SetTitle("AdcPerSigLin"); hAdcPerSigLinReCal1->GetXaxis()->SetTitle("AdcPerSigLin"); hAdcPerSigLinReCal1->GetXaxis()->CenterTitle(); //hAdcPerSigLinReCal1->SetBit(TH1::kCanRebin); hSigLinPerSigCorReCal1=new TH1F("hSigLinPerSigCorReCal1", "hSigLinPerSigCorReCal1", 1100,-1,10); hSigLinPerSigCorReCal1->SetFillColor(0); hSigLinPerSigCorReCal1->SetTitle("SigLinPerSigCor"); hSigLinPerSigCorReCal1->GetXaxis()->SetTitle("SigLinPerSigCor"); hSigLinPerSigCorReCal1->GetXaxis()->CenterTitle(); //hSigLinPerSigCorReCal1->SetBit(TH1::kCanRebin); hSigCorPerSigMapReCal1=new TH1F("hSigCorPerSigMapReCal1", "hSigCorPerSigMapReCal1", 10000,-1,100); hSigCorPerSigMapReCal1->SetFillColor(0); hSigCorPerSigMapReCal1->SetTitle("SigCorPerSigMap"); hSigCorPerSigMapReCal1->GetXaxis()->SetTitle("SigCorPerSigMap"); hSigCorPerSigMapReCal1->GetXaxis()->CenterTitle(); //hSigCorPerSigMapReCal1->SetBit(TH1::kCanRebin); hSigMapPerMipReCal1=new TH1F("hSigMapPerMipReCal1", "hSigMapPerMipReCal1", 10100,-10,1000); hSigMapPerMipReCal1->SetFillColor(0); hSigMapPerMipReCal1->SetTitle("SigMapPerMip"); hSigMapPerMipReCal1->GetXaxis()->SetTitle("SigMapPerMip"); hSigMapPerMipReCal1->GetXaxis()->CenterTitle(); //hSigMapPerMipReCal1->SetBit(TH1::kCanRebin); hMipPerGeVReCal1=new TH1F("hMipPerGeVReCal1","hMipPerGeVReCal1", 10100,-10,1000); hMipPerGeVReCal1->SetFillColor(0); hMipPerGeVReCal1->SetTitle("MipPerGeV"); hMipPerGeVReCal1->GetXaxis()->SetTitle("MipPerGeV"); hMipPerGeVReCal1->GetXaxis()->CenterTitle(); //hMipPerGeVReCal1->SetBit(TH1::kCanRebin); hAdcPerPeReCal1=new TH1F("hAdcPerPeReCal1","hAdcPerPeReCal1", 8080,-20,2000); hAdcPerPeReCal1->SetFillColor(0); hAdcPerPeReCal1->SetTitle("AdcPerPe"); hAdcPerPeReCal1->GetXaxis()->SetTitle("AdcPerPe"); hAdcPerPeReCal1->GetXaxis()->CenterTitle(); //hAdcPerPeReCal1->SetBit(TH1::kCanRebin); hAdcPerSigDrfReCal1=new TH1F("hAdcPerSigDrfReCal1", "hAdcPerSigDrfReCal1", 1100,-1,10); hAdcPerSigDrfReCal1->SetFillColor(0); hAdcPerSigDrfReCal1->SetTitle("AdcPerSigDrf"); hAdcPerSigDrfReCal1->GetXaxis()->SetTitle("AdcPerSigDrf"); hAdcPerSigDrfReCal1->GetXaxis()->CenterTitle(); //hAdcPerSigDrfReCal1->SetBit(TH1::kCanRebin); hAdcPerSigLinOnlyReCal1=new TH1F("hAdcPerSigLinOnlyReCal1", "hAdcPerSigLinOnlyReCal1", 2100,-1,20); hAdcPerSigLinOnlyReCal1->SetFillColor(0); hAdcPerSigLinOnlyReCal1->SetTitle("AdcPerSigLinOnly"); hAdcPerSigLinOnlyReCal1->GetXaxis()->SetTitle("AdcPerSigLinOnly"); hAdcPerSigLinOnlyReCal1->GetXaxis()->CenterTitle(); //hAdcPerSigLinOnlyReCal1->SetBit(TH1::kCanRebin); //stripend 2 hAdcPerSigLinReCal2=new TH1F("hAdcPerSigLinReCal2","hAdcPerSigLinReCal2", 1100,-1,10); hAdcPerSigLinReCal2->SetFillColor(0); hAdcPerSigLinReCal2->SetTitle("AdcPerSigLin"); hAdcPerSigLinReCal2->GetXaxis()->SetTitle("AdcPerSigLin"); hAdcPerSigLinReCal2->GetXaxis()->CenterTitle(); //hAdcPerSigLinReCal2->SetBit(TH1::kCanRebin); hSigLinPerSigCorReCal2=new TH1F("hSigLinPerSigCorReCal2", "hSigLinPerSigCorReCal2", 1100,-1,10); hSigLinPerSigCorReCal2->SetFillColor(0); hSigLinPerSigCorReCal2->SetTitle("SigLinPerSigCor"); hSigLinPerSigCorReCal2->GetXaxis()->SetTitle("SigLinPerSigCor"); hSigLinPerSigCorReCal2->GetXaxis()->CenterTitle(); //hSigLinPerSigCorReCal2->SetBit(TH1::kCanRebin); hSigCorPerSigMapReCal2=new TH1F("hSigCorPerSigMapReCal2", "hSigCorPerSigMapReCal2", 10000,-1,100); hSigCorPerSigMapReCal2->SetFillColor(0); hSigCorPerSigMapReCal2->SetTitle("SigCorPerSigMap"); hSigCorPerSigMapReCal2->GetXaxis()->SetTitle("SigCorPerSigMap"); hSigCorPerSigMapReCal2->GetXaxis()->CenterTitle(); //hSigCorPerSigMapReCal2->SetBit(TH1::kCanRebin); hSigMapPerMipReCal2=new TH1F("hSigMapPerMipReCal2", "hSigMapPerMipReCal2", 10100,-10,1000); hSigMapPerMipReCal2->SetFillColor(0); hSigMapPerMipReCal2->SetTitle("SigMapPerMip"); hSigMapPerMipReCal2->GetXaxis()->SetTitle("SigMapPerMip"); hSigMapPerMipReCal2->GetXaxis()->CenterTitle(); //hSigMapPerMipReCal2->SetBit(TH1::kCanRebin); hMipPerGeVReCal2=new TH1F("hMipPerGeVReCal2","hMipPerGeVReCal2", 10100,-10,1000); hMipPerGeVReCal2->SetFillColor(0); hMipPerGeVReCal2->SetTitle("MipPerGeV"); hMipPerGeVReCal2->GetXaxis()->SetTitle("MipPerGeV"); hMipPerGeVReCal2->GetXaxis()->CenterTitle(); //hMipPerGeVReCal2->SetBit(TH1::kCanRebin); hAdcPerPeReCal2=new TH1F("hAdcPerPeReCal2","hAdcPerPeReCal2", 8080,-20,2000); hAdcPerPeReCal2->SetFillColor(0); hAdcPerPeReCal2->SetTitle("AdcPerPe"); hAdcPerPeReCal2->GetXaxis()->SetTitle("AdcPerPe"); hAdcPerPeReCal2->GetXaxis()->CenterTitle(); //hAdcPerPeReCal2->SetBit(TH1::kCanRebin); hAdcPerSigDrfReCal2=new TH1F("hAdcPerSigDrfReCal2", "hAdcPerSigDrfReCal2", 1100,-1,10); hAdcPerSigDrfReCal2->SetFillColor(0); hAdcPerSigDrfReCal2->SetTitle("AdcPerSigLinDrf"); hAdcPerSigDrfReCal2->GetXaxis()->SetTitle("AdcPerSigLinDrf"); hAdcPerSigDrfReCal2->GetXaxis()->CenterTitle(); //hAdcPerSigDrfReCal2->SetBit(TH1::kCanRebin); hAdcPerSigLinOnlyReCal2=new TH1F("hAdcPerSigLinOnlyReCal2", "hAdcPerSigLinOnlyReCal2", 2100,-1,20); hAdcPerSigLinOnlyReCal2->SetFillColor(0); hAdcPerSigLinOnlyReCal2->SetTitle("AdcPerSigLinOnly"); hAdcPerSigLinOnlyReCal2->GetXaxis()->SetTitle("AdcPerSigLinOnly"); hAdcPerSigLinOnlyReCal2->GetXaxis()->CenterTitle(); //hAdcPerSigLinOnlyReCal2->SetBit(TH1::kCanRebin); } Calibrator& cal=Calibrator::Instance(); //VldTimeStamp (const time_t &t, const Int_t nsec) //VldTimeStamp vldts; VldTimeStamp vldts(s.TimeSec,0); VldContext vc(static_cast<Detector::Detector_t>(s.Detector), static_cast<SimFlag::ESimFlag>(s.SimFlag),vldts); cal.ReInitialise(vc);//same as reset //vc.Print(); //cal.PrintConfig(cout); MAXMSG("MeuAnalysis",logLevel,1000) <<"Recalibrate: track "<<s.VtxPlane<<"->"<<s.EndPlane<<endl; //zero all the current values for (map<Int_t,MeuHitInfo>::iterator it=plInfo.begin(); it!=plInfo.end();++it){ //cache the current cp MeuHitInfo& cp=it->second; MAXMSG("MeuAnalysis",logLevel,1000) <<"pl="<<cp.Plane <<", st="<<cp.Strip <<", sigCor="<<cp.SigCor <<", sigMap="<<cp.SigMap <<", lpos="<<cp.LPos<<endl; cp.SigMap=0; cp.SigMap1=0; cp.SigMap2=0; cp.MCSigMap=0; if (fRecalibrateSigCor) { MAXMSG("MeuAnalysis",Msg::kInfo,1) <<"Reapplying strip-to-strip calibration" <<", s.SimFlag="<<s.SimFlag<<endl; cp.SigCor=0; cp.SigCor1=0; cp.SigCor2=0; cp.SigCorTrk1=0; cp.SigCorTrk2=0; } else { MAXMSG("MeuAnalysis",Msg::kInfo,1) <<"NOT Reapplying strip-to-strip calibration"<<endl; } if (fRecalibrateSigLin) { MAXMSG("MeuAnalysis",Msg::kInfo,1) <<"Reapplying SigLin calibration"<<endl; cp.SigLin=0; cp.SigLin1=0; cp.SigLin2=0; } else { MAXMSG("MeuAnalysis",Msg::kInfo,1) <<"NOT Reapplying SigLin calibration (hence no drift)"<<endl; } if (fRecalibratePE) { MAXMSG("MeuAnalysis",Msg::kInfo,1) <<"Reapplying PE calibration"<<endl; cp.Pe=0; cp.Pe1=0; cp.Pe2=0; } else { MAXMSG("MeuAnalysis",Msg::kInfo,1) <<"NOT Reapplying PE calibration"<<endl; } //these should already be zero... but to be on the safe side: cp.SigLinOnly=0; cp.SigLinOnly1=0; cp.SigLinOnly2=0; cp.SigDrf=0; cp.SigDrf1=0; cp.SigDrf2=0; } TClonesArray& slcTca=(*ntp.slc); TClonesArray& stpTca=(*ntp.stp); //get the slice associated with this event const NtpSRSlice* pslc= dynamic_cast<NtpSRSlice*>(slcTca[evt.slc]); const NtpSRSlice& slc=(*pslc); //loop over strips in slc for (Int_t i=0;i<slc.nstrip;++i) { const NtpSRStrip* pstp= dynamic_cast<NtpSRStrip*>(stpTca[slc.stp[i]]); const NtpSRStrip& stp=(*pstp); Int_t pl=stp.plane; //check if the plane is in the track //this cut is not in NuCuts::ExtractPlInfo but it should //never matter since the track window is always within the track if (!this->PlaneIsWithinTrack(s,pl, s.VtxPlane,s.EndPlane)) continue; //filter out the times that are way off to make the //Near and Far detectors more similar //get the time0, if between 0 and 1 switch to time1 Double_t time=stp.time0; if (time<0 || time>1) time=stp.time1; if (time<0 || time>1) cout<<"Ahhhhh, time"<<endl;; if (time>s.MedianTime+(400*Munits::ns) || time<s.MedianTime-(400*Munits::ns)) { MAXMSG("MeuCuts",Msg::kInfo,10) <<"Cut strip: time="<<time <<", pl="<<stp.plane<<", st="<<stp.strip <<", sigCor="<<stp.ph0.sigcor+stp.ph1.sigcor <<", dT="<<(time-s.MedianTime)/Munits::ns<<endl; continue;//messes up beg/end planes (if not filtered bad trks)! } //cache the current cp MeuHitInfo& cp=plInfo[pl]; //get the raw adcs Float_t adc1=stp.ph0.raw;//east Float_t adc2=stp.ph1.raw;//west PlexStripEndId seidE(static_cast<Detector::Detector_t> (s.Detector),pl,stp.strip, StripEnd::kEast); PlexStripEndId seidW(static_cast<Detector::Detector_t> (s.Detector),pl,stp.strip, StripEnd::kWest); //sigDrf and sigLinOnly values are NOT used to get the MEU number //get the sigLinOnly value (applied to strip) //useful for systematic error evaluation //NOTE: in macro have to set: //cal.LinearityCalScheme().Set("BucketMode=0"); //for the results to have any meaning Float_t sigLinOnly1=0; Float_t sigLinOnly2=0; if (adc1>0) sigLinOnly1=cal.GetLinearized(adc1,seidE); else sigLinOnly1=0; if (adc2>0) sigLinOnly2=cal.GetLinearized(adc2,seidW); else sigLinOnly2=0; //now sum up everything in the plane cp.SigLinOnly+=sigLinOnly1+sigLinOnly2; cp.SigLinOnly1+=sigLinOnly1; cp.SigLinOnly2+=sigLinOnly2; //get the drift correction applied to adcs Float_t sigDrf1=0; Float_t sigDrf2=0; if (adc1>0) sigDrf1=cal.GetDriftCorrected(adc1,seidE); else sigDrf1=0; if (adc2>0) sigDrf2=cal.GetDriftCorrected(adc2,seidW); else sigDrf2=0; //now sum up everything in the plane cp.SigDrf+=sigDrf1+sigDrf2; cp.SigDrf1+=sigDrf1; cp.SigDrf2+=sigDrf2; //section to recalibrate PEs if required Float_t pe1=stp.ph0.pe;//east Float_t pe2=stp.ph1.pe;//west if (fRecalibratePE){ if (adc1>0) pe1=cal.GetPhotoElectrons(adc1,seidE); else pe1=0; if (adc2>0) pe2=cal.GetPhotoElectrons(adc2,seidW); else pe2=0; //now sum up everything in the plane cp.Pe+=pe1+pe2; cp.Pe1+=pe1; cp.Pe2+=pe2; } //section to recalibrate siglin if required Float_t sigLin1=stp.ph0.siglin;//east Float_t sigLin2=stp.ph1.siglin;//west if (fRecalibrateSigLin){ //linearity calibration is done on digits (not strips) //this matters for the ND! //It is not totally correct to recalibrate with the standard //ntuples because the digit information is not available for ND //however the calibrator can implement an accurate fudge: //cal.LinearityCalScheme().Set("BucketMode=0"); MAXMSG("MeuAnalysis",Msg::kInfo,1) <<"Recalibrating drift and linearity for detector=" <<s.Detector<<endl; //first apply the linearity if (adc1>0) sigLin1=cal.GetLinearized(adc1,seidE); else sigLin1=0; if (adc2>0) sigLin2=cal.GetLinearized(adc2,seidW); else sigLin2=0; //then do the drift correction (to sigLin - driftCor this time) if (sigLin1>0) sigLin1=cal.GetDriftCorrected(sigLin1,seidE); else sigLin1=0; if (sigLin2>0) sigLin2=cal.GetDriftCorrected(sigLin2,seidW); else sigLin2=0; //now sum up everything in the plane cp.SigLin+=sigLin1+sigLin2; cp.SigLin1+=sigLin1; cp.SigLin2+=sigLin2; } //section to recalibrate sigcor if required Float_t sigCor1=stp.ph0.sigcor;//east Float_t sigCor2=stp.ph1.sigcor;//west if (fRecalibrateSigCor){ MAXMSG("MeuAnalysis",Msg::kInfo,1) <<"Recalibrating strip-to-strip..."<<endl; if (s.Detector!=Detector::kNear) { if (sigLin1>0) sigCor1=cal.GetStripToStripCorrected(sigLin1, seidE); } //else sigCor1=0; if (sigLin2>0) sigCor2=cal.GetStripToStripCorrected(sigLin2, seidW); else sigCor2=0; //sum both stripends Float_t sigCorBoth=sigCor1+sigCor2; //sum the sigcor in the object cp.SigCor+=sigCorBoth; cp.SigCor1+=sigCor1; cp.SigCor2+=sigCor2; //store the sigcor of both ends of the tracked strip //this means that fSigCor!=fSigCor1+fSigCor2 always if (stp.strip==cp.Strip){ cp.SigCorTrk1=sigCor1; cp.SigCorTrk2=sigCor2; } } else if (fRecalibrateSigLin) { MAXMSG("MeuAnalysis",Msg::kInfo,1) <<"Using ntp strip-to-strip calibration on recalibrated" <<" SigLin values..."<<endl; //if siglin was recalibrated then have to feed the value through //the value of the sigcor correction in the ntp is just //stp.ph0.sigcor/stp.ph0.siglin so multiply by that: sigCor1=sigLin1*stp.ph0.sigcor/stp.ph0.siglin;//east sigCor2=sigLin2*stp.ph1.sigcor/stp.ph1.siglin;//west } //find the lpos for this plane then calibrate Float_t lpos=cp.LPos; Float_t sigMap1=0; Float_t sigMap2=0; //only calibrate the east side for !ND //also only do the calibration if there is a hit on stripend if (s.Detector!=Detector::kNear) { if (sigCor1>0) sigMap1=cal.GetAttenCorrectedTpos(sigCor1, lpos,seidE); } if (sigCor2>0) sigMap2=cal.GetAttenCorrectedTpos(sigCor2, lpos,seidW); Float_t sigMapBoth=sigMap1+sigMap2; //sum the sigmap cp.SigMap+=sigMapBoth; cp.SigMap1+=sigMap1; cp.SigMap2+=sigMap2; MAXMSG("MeuAnalysis",logLevel,1000) <<"strip="<<stp.strip<<", pl="<<pl<<", lpos="<<lpos<<endl <<" sigCor1="<<sigCor1<<", sigMap1="<<sigMap1<<endl <<" sigCor2="<<sigCor2<<", sigMap2="<<sigMap2<<endl <<endl; //calibrate "true MC" values as well if (s.SimFlag==SimFlag::kMC){ //get the MC lpos Float_t lpos=cp.MCLPos; Float_t sigMap1=0; Float_t sigMap2=0; if (sigCor1>0) { sigMap1=cal.GetAttenCorrectedTpos(sigCor1,lpos,seidE); } if (sigCor2>0) { sigMap2=cal.GetAttenCorrectedTpos(sigCor2,lpos,seidW); } //calc sigmap and add to cp sigMapBoth=sigMap1+sigMap2; cp.MCSigMap+=sigMapBoth; } //fill calibration value histos //can't do mip or gev ones if (sigLin1) hAdcPerSigLinReCal1->Fill(adc1/sigLin1); if (sigCor1) hSigLinPerSigCorReCal1->Fill(sigLin1/sigCor1); if (sigCor1) hSigCorPerSigMapReCal1->Fill(sigCor1/sigMap1); if (pe1) hAdcPerPeReCal1->Fill(adc1/pe1); if (sigLinOnly1) hAdcPerSigLinOnlyReCal1->Fill(adc1/sigLinOnly1); if (sigDrf1) hAdcPerSigDrfReCal1->Fill(adc1/sigDrf1); //stripend 2 if (sigLin2) hAdcPerSigLinReCal2->Fill(adc2/sigLin2); if (sigCor2) hSigLinPerSigCorReCal2->Fill(sigLin2/sigCor2); if (sigCor2) hSigCorPerSigMapReCal2->Fill(sigCor2/sigMap2); if (pe2) hAdcPerPeReCal2->Fill(adc2/pe2); if (sigLinOnly2) hAdcPerSigLinOnlyReCal2->Fill(adc2/sigLinOnly2); if (sigDrf2) hAdcPerSigDrfReCal2->Fill(adc2/sigDrf2); //print out static Bool_t print=false; if (print) { if (sigDrf2) { if (adc2/sigDrf2>2) { MAXMSG("MeuAnalysis",Msg::kInfo,1000) <<"adc2/sigDrf2="<<adc2/sigDrf2 <<", adc2="<<adc2<<", sigDrf2="<<sigDrf2<<endl; } } if (sigLinOnly2) { if (adc2/sigLinOnly2>2) { MAXMSG("MeuAnalysis",Msg::kInfo,1000) <<"adc2/sigLinOnly2="<<adc2/sigLinOnly2 <<", adc2="<<adc2<<", sigLinOnly2="<<sigLinOnly2 <<", pl,st="<<cp.Plane<<","<<cp.Strip<<endl; } } } }//loop over strips in slc //print Bool_t printMap=false; if (printMap){ for (map<Int_t,MeuHitInfo>::iterator it=plInfo.begin(); it!=plInfo.end();++it){ //cache the current cp MeuHitInfo& cp=it->second; MAXMSG("MeuAnalysis",Msg::kInfo,1000) <<"Printing: pl="<<cp.Plane <<", st="<<cp.Strip <<", sigCor="<<cp.SigCor <<", sigMap="<<cp.SigMap <<", lpos="<<cp.LPos<<endl; } } }

void MeuAnalysis::RecalibrateAtNu (  AtmosEvent &  ev, MeuSummary &  s, std::map< Int_t, MeuHitInfo > &  plInfo )  const [private]

Definition at line 2666 of file MeuAnalysis.cxx. References Calibrator::GetAttenCorrectedTpos(), Calibrator::GetDriftCorrected(), Calibrator::GetLinearized(), Calibrator::GetPhotoElectrons(), Calibrator::GetStripToStripCorrected(), Calibrator::Instance(), MeuHitInfo::LPos, MAXMSG, MeuHitInfo::MCLPos, MeuHitInfo::MCSigMap, MeuHitInfo::Pe, MeuHitInfo::Pe1, MeuHitInfo::Pe2, AtmosStrip::Plane, MeuHitInfo::Plane, Calibrator::PrintConfig(), AtmosStrip::Qadc, AtmosStrip::QPE, Calibrator::ReInitialise(), s(), CfgPromptConfigurable::Set(), MeuHitInfo::SigCor, MeuHitInfo::SigCor1, MeuHitInfo::SigCor2, MeuHitInfo::SigCorTrk1, MeuHitInfo::SigCorTrk2, MeuHitInfo::SigLin, MeuHitInfo::SigLin1, MeuHitInfo::SigLin2, MeuHitInfo::SigMap, MeuHitInfo::SigMap1, MeuHitInfo::SigMap2, AtmosStrip::Strip, and MeuHitInfo::Strip. Referenced by MakeSummaryTreeWithAtNu(). { //variable to turn on all the useful messages if required Msg::LogLevel_t logLevel=Msg::kDebug; //Msg::LogLevel_t logLevelSts=Msg::kDebug; //Msg::LogLevel_t logLevelLin=Msg::kInfo; Bool_t fRecalibrateSigLin=true;//have to configure the calibrator too! Bool_t fRecalibrateSigCor=true;//have to configure the calibrator too! Bool_t fRecalibratePE=true;//have to configure the calibrator too! Calibrator& cal=Calibrator::Instance(); //VldTimeStamp (const time_t &t, const Int_t nsec) //VldTimeStamp vldts; VldTimeStamp vldts(ev.UnixTime,0); VldContext vc(static_cast<Detector::Detector_t>(s.Detector), static_cast<SimFlag::ESimFlag>(s.SimFlag),vldts); cal.ReInitialise(vc); static Bool_t firstTime=false;//switch this off with false if (firstTime){ if (s.Detector==Detector::kFar){ cal.Set("LinCalibrator=PulserLinearityCalScheme"); } else if (s.Detector==Detector::kNear){ cal.Set("LinCalibrator=QuadLinearityCalScheme"); } else cout<<"Ahh, detector="<<s.Detector<<endl; cout<<"Actually using calibrator configured as:"<<endl; cal.PrintConfig(cout); firstTime=false; } MAXMSG("MeuAnalysis",logLevel,1000) <<"Recalibrate: track "<<s.VtxPlane<<"->"<<s.EndPlane<<endl; //zero all the current values for (map<Int_t,MeuHitInfo>::iterator it=plInfo.begin(); it!=plInfo.end();++it){ //cache the current cp MeuHitInfo& cp=it->second; MAXMSG("MeuAnalysis",logLevel,1000) <<"pl="<<cp.Plane <<", st="<<cp.Strip <<", sigCor="<<cp.SigCor <<", sigMap="<<cp.SigMap <<", lpos="<<cp.LPos<<endl; cp.SigMap=0; cp.SigMap1=0; cp.SigMap2=0; cp.MCSigMap=0; if (fRecalibrateSigCor) { MAXMSG("MeuAnalysis",Msg::kInfo,1) <<"Reapplying strip-to-strip calibration, cp.SimFlag=" <<s.SimFlag<<endl; cp.SigCor=0; cp.SigCor1=0; cp.SigCor2=0; cp.SigCorTrk1=0; cp.SigCorTrk2=0; } else { MAXMSG("MeuAnalysis",Msg::kInfo,1) <<"NOT Reapplying strip-to-strip calibration"<<endl; } if (fRecalibrateSigLin) { MAXMSG("MeuAnalysis",Msg::kInfo,1) <<"Reapplying SigLin calibration"<<endl; cp.SigLin=0; cp.SigLin1=0; cp.SigLin2=0; } else { MAXMSG("MeuAnalysis",Msg::kInfo,1) <<"NOT Reapplying SigLin calibration"<<endl; } if (fRecalibratePE) { MAXMSG("MeuAnalysis",Msg::kInfo,1) <<"Reapplying PE calibration"<<endl; cp.Pe=0; cp.Pe1=0; cp.Pe2=0; } else { MAXMSG("MeuAnalysis",Msg::kInfo,1) <<"NOT Reapplying PE calibration"<<endl; } } TClonesArray& strips=(*ev.StripList); Int_t numStrips=strips.GetEntries(); //loop over the strips in the snarl for (Int_t hit=0;hit<numStrips;hit++){ const AtmosStrip* strip=dynamic_cast<AtmosStrip*>(strips[hit]); const AtmosStrip& stp=(*strip); Int_t pl=stp.Plane; //check if the plane is in the track if (!this->PlaneIsWithinTrack(s,pl, s.VtxPlane,s.EndPlane)) continue; //cache the current cp MeuHitInfo& cp=plInfo[pl]; Float_t adc1=stp.Qadc[0];//east Float_t adc2=stp.Qadc[1];//west Float_t pe1=stp.QPE[0];//east Float_t pe2=stp.QPE[1];//west Float_t sigLin1=stp.Qadc[0];//east//siglin is no longer in ntp Float_t sigLin2=stp.Qadc[1];//west//so use adcs instead Float_t sigCor1=stp.Qadc[0];//east//siglin is no longer in ntp Float_t sigCor2=stp.Qadc[1];//west//so use adcs instead //Float_t sigCor1=stp.QPEcorr[0];//east//this is pe size so don't use it! //Float_t sigCor2=stp.QPEcorr[1];//west PlexStripEndId seidE(static_cast<Detector::Detector_t> (s.Detector),pl,stp.Strip, StripEnd::kEast); PlexStripEndId seidW(static_cast<Detector::Detector_t> (s.Detector),pl,stp.Strip, StripEnd::kWest); //section to recalibrate PEs if required if (fRecalibratePE){ if (adc1>0) pe1=cal.GetPhotoElectrons(adc1,seidE); else pe1=0; if (adc2>0) pe2=cal.GetPhotoElectrons(adc2,seidW); else pe2=0; //now sum up everything in the plane cp.Pe+=pe1+pe2; cp.Pe1+=pe1; cp.Pe2+=pe2; } //section to recalibrate siglin if required if (fRecalibrateSigLin){ //first do the drift correction if (adc1>0) sigLin1=cal.GetDriftCorrected(adc1,seidE); else sigLin1=0; if (adc2>0) sigLin2=cal.GetDriftCorrected(adc2,seidW); else sigLin2=0; //then apply the linearity (to sigLin - driftCor this time) if (adc1>0) sigLin1=cal.GetLinearized(sigLin1,seidE); else sigLin1=0; if (adc2>0) sigLin2=cal.GetLinearized(sigLin2,seidW); else sigLin2=0; //now sum up everything in the plane cp.SigLin+=sigLin1+sigLin2; cp.SigLin1+=sigLin1; cp.SigLin2+=sigLin2; } //section to recalibrate sigcor if required // have to do this for AtNu for MC AND DATA // since it has pe-sized units of "cor" if (fRecalibrateSigCor){// && s.SimFlag==SimFlag::kData){ if (sigLin1>0) sigCor1=cal.GetStripToStripCorrected(sigLin1, seidE); else sigCor1=0; if (sigLin2>0) sigCor2=cal.GetStripToStripCorrected(sigLin2, seidW); else sigCor2=0; //sum the sigcor in the object cp.SigCor+=sigCor1+sigCor2; cp.SigCor1+=sigCor1; cp.SigCor2+=sigCor2; //store the sigcor of both ends of the tracked strip //this means that fSigCor!=fSigCor1+fSigCor2 always if (stp.Strip==cp.Strip){ cp.SigCorTrk1=sigCor1; cp.SigCorTrk2=sigCor2; } } //find the lpos for this plane then calibrate Float_t lpos=cp.LPos; Float_t sigMap1=0; Float_t sigMap2=0; //only do the look up if non zero if (sigCor1>0) sigMap1=cal.GetAttenCorrectedTpos(sigCor1, lpos,seidE); if (sigCor2>0) sigMap2=cal.GetAttenCorrectedTpos(sigCor2, lpos,seidW); Float_t sigMapBoth=sigMap1+sigMap2; //sum the sigmap cp.SigMap+=sigMapBoth; //store the sigmap of both ends of the tracked strip //this means that fSigMap!=fSigMap1+fSigMap2 always //if (st.Strip==cp.Strip){ cp.SigMap1+=sigMap1; cp.SigMap2+=sigMap2; //} MAXMSG("MeuAnalysis",logLevel,1000) <<"strip="<<stp.Strip<<", pl="<<pl<<", lpos="<<lpos<<endl <<" sigCor1="<<sigCor1<<", sigMap1="<<sigMap1<<endl <<" sigCor2="<<sigCor2<<", sigMap2="<<sigMap2<<endl <<endl; //do MC as well if (s.SimFlag==SimFlag::kMC){ //get the MC lpos Float_t lpos=cp.MCLPos; Float_t sigMap1=0; Float_t sigMap2=0; if (sigCor1>0) { sigMap1=cal.GetAttenCorrectedTpos(sigCor1,lpos,seidE); } if (sigCor2>0) { sigMap2=cal.GetAttenCorrectedTpos(sigCor2,lpos,seidW); } //calc sigmap and add to cp sigMapBoth=sigMap1+sigMap2; cp.MCSigMap+=sigMapBoth; } } //print Bool_t printMap=false; if (printMap){ for (map<Int_t,MeuHitInfo>::iterator it=plInfo.begin(); it!=plInfo.end();++it){ //cache the current cp MeuHitInfo& cp=it->second; MAXMSG("MeuAnalysis",Msg::kInfo,1000) <<"Printing: pl="<<cp.Plane <<", st="<<cp.Strip <<", sigCor="<<cp.SigCor <<", sigMap="<<cp.SigMap <<", lpos="<<cp.LPos<<endl; } } }

void MeuAnalysis::RecalibratePE (  Bool_t  recalibrate = true  )  [static]

Definition at line 201 of file MeuAnalysis.cxx. References fRecalibratePE, and MSG. { MSG("MeuAnalysis",Msg::kInfo) <<"Setting fRecalibratePE="<<recalibrate<<endl; fRecalibratePE=recalibrate; if (recalibrate) { MSG("MeuAnalysis",Msg::kInfo) <<"Note that the calibrator has to be correctly configured" <<" for the recalibration to work!"<<endl; } }

void MeuAnalysis::RecalibrateSigCor (  Bool_t  recalibrate = true  )  [static]

Definition at line 187 of file MeuAnalysis.cxx. References fRecalibrateSigCor, and MSG. { MSG("MeuAnalysis",Msg::kInfo) <<"Setting fRecalibrateSigCor="<<recalibrate<<endl; fRecalibrateSigCor=recalibrate; if (recalibrate) { MSG("MeuAnalysis",Msg::kInfo) <<"Note that the calibrator has to be correctly configured" <<" for the recalibration to work!"<<endl; } }

void MeuAnalysis::RecalibrateSigLin (  Bool_t  recalibrate = true  )  [static]

Definition at line 173 of file MeuAnalysis.cxx. References fRecalibrateSigLin, and MSG. { MSG("MeuAnalysis",Msg::kInfo) <<"Setting fRecalibrateSigLin="<<recalibrate<<endl; fRecalibrateSigLin=recalibrate; if (recalibrate) { MSG("MeuAnalysis",Msg::kInfo) <<"Note that the calibrator has to be correctly configured" <<" for the recalibration to work!"<<endl; } }

void MeuAnalysis::SetLoopVariables (  Int_t  event, Int_t  printMode = 1 )  [private]

Definition at line 743 of file MeuAnalysis.cxx. References fChain, fChainBD, fEntries, and MSG. Referenced by BasicPlots(), BasicReco(), MakeSummaryTreeWithAtNu(), MakeSummaryTreeWithNtpSt(), N_1Plots(), SnarlList(), SpillPlots(), and TestEventLoop(). { if (printMode==1){ Float_t fract=ceil(fEntries/20.); if (ceil(((Float_t)entry)/fract)==((Float_t)entry)/fract){ MSG("MeuAnalysis",Msg::kInfo) <<"Fraction of loop complete: "<<entry <<"/"<<fEntries<<" (" <<(Int_t)(100.*entry/fEntries)<<"%)"<<endl; } } //get the snarl/entry fChain->GetEvent(entry); if (fChainBD) fChainBD->GetEvent(entry); }

void MeuAnalysis::SetRootFileObjects (   )  [private]

Definition at line 593 of file MeuAnalysis.cxx. References NtpStRecord::evthdr, fChain, fChainBD, fEntries, fEntriesBD, fev, fRec, fRecBD, RecCandHeader::GetEvent(), MSG, NtpSREventSummary::nslice, and NtpBDLiteRecord::tortgt. { MSG("MeuAnalysis",Msg::kInfo) <<"Running the SetRootFileObjects method..."<<endl; if (fUseAtNu) { fChain->SetBranchAddress("evt", &fev); //set the number of events in the tree fEntries=static_cast<Int_t>(fChain->GetEntries()); MSG("MeuAnalysis",Msg::kInfo) <<"AtNuTree has "<<fEntries<<" entries"<<endl; if (fEntries>0){ //get the snarl/event fChain->GetEvent(0); //fRunNumber=(*fev).Run; MSG("MeuAnalysis",Msg::kInfo) <<"SetRootFileObjects: first UnixTime="<<(*fev).UnixTime<<endl; } } else { //set up tree fChain->SetBranchAddress("NtpStRecord",&fRec); if (fChainBD) fChainBD->SetBranchAddress("NtpBDLiteRecord",&fRecBD); //fBDtortgt //set the number of events in the tree fEntries=static_cast<Int_t>(fChain->GetEntries()); if (fChainBD) fEntriesBD=static_cast<Int_t>(fChainBD->GetEntries()); MSG("MeuAnalysis",Msg::kInfo) <<"NtpSt tree has "<<fEntries<<" entries"<<endl; MSG("MeuAnalysis",Msg::kInfo) <<"NtpBDLite tree has "<<fEntriesBD<<" entries"<<endl; if (fEntries>0){ //get the snarl/event fChain->GetEvent(0); //MSG("MeuAnalysis",Msg::kInfo) //<<"First run number="<<fRec->fHeader.fRun<<endl; MSG("MeuAnalysis",Msg::kInfo) <<"SetRootFileObjects: first snarl nslice=" <<fRec->evthdr.nslice<<endl; } if (fEntriesBD>0){ //get the snarl/event fChainBD->GetEvent(0); MSG("MeuAnalysis",Msg::kInfo) <<"First tortgt="<<fRecBD->tortgt<<" E12 POT"<<endl; } } MSG("MeuAnalysis",Msg::kDebug) <<"Finished the SetRootFileObjects method"<<endl; }

void MeuAnalysis::SnarlList (   )

Definition at line 3183 of file MeuAnalysis.cxx. References abs(), MeuReco::CalcFCPC(), MeuReco::CalcVtx(), MeuReco::CalcVtxSpecialND(), MeuSummary::DistToEdgeFid, NtpSRTrack::end, MeuSummary::EndPlane, MeuSummary::Event, NtpStRecord::evt, MeuCuts::ExtractPlInfo(), MeuCuts::FillSTSumDetails(), RecRecordImp< T >::GetHeader(), RecDataHeader::GetRun(), RecPhysicsHeader::GetSnarl(), InitialiseLoopVariables(), MAXMSG, MSG, OpenTxtFile(), MeuSummary::PC, NtpSRVertex::plane, MeuSummary::Reset(), MeuSummary::RFid, MeuSummary::Run, s(), SetLoopVariables(), MeuSummary::Snarl, NtpStRecord::trk, and NtpSRTrack::vtx. { MSG("MeuAnalysis",Msg::kInfo) <<" ** Running SnarlList method... **"<<endl; ofstream& snarlList=*(this->OpenTxtFile(100,"snarlList")); //reconstruction object MeuReco reco; MeuCuts cuts; MeuSummary s; //variable to turn on all the useful messages if required Msg::LogLevel_t logLevel=Msg::kDebug; this->InitialiseLoopVariables(); //for(Int_t entry=55000;entry<fEntries;entry++){ for(Int_t entry=0;entry<fEntries;entry++){ this->SetLoopVariables(entry); MSG("MeuAnalysis",logLevel) <<"Entry="<<entry<<endl; //make a reference instead of a pointer NtpStRecord& ntp=(*fRec); const RecCandHeader& rec=ntp.GetHeader(); MAXMSG("MeuAnalysis",Msg::kInfo,1) <<"First entry: run="<<rec.GetRun() <<", snarl="<<rec.GetSnarl()<<endl; //MeuSummary s; s.Reset(); s.RFid=1.0; s.DistToEdgeFid=0.4; s.Event=entry; cuts.FillSTSumDetails(ntp,s,-1,-1); //a map to hold the info for each plane in the entry map<Int_t,MeuHitInfo> plInfo; //ensure there is only 1 track in the entry TClonesArray& trkTca=(*ntp.trk); Int_t numTrks=trkTca.GetEntries(); if (numTrks!=1) continue; const NtpSRTrack* ptrk=dynamic_cast<NtpSRTrack*>(trkTca[0]); const NtpSRTrack& trk=(*ptrk); //cut on number of planes Int_t numPlanes=abs(trk.vtx.plane-trk.end.plane); if (numPlanes<6) continue; TClonesArray& evtTca=(*ntp.evt); const NtpSREvent* pevt=dynamic_cast<NtpSREvent*>(evtTca[0]); cuts.ExtractPlInfo(ntp,s,plInfo,*pevt); reco.CalcVtx(s,plInfo); //check if track stops in spectrometer if (s.EndPlane>120) continue; reco.CalcFCPC(s,plInfo); //recalculate the vtx of the tracks coming in from the spectrometer reco.CalcVtxSpecialND(s,plInfo); if (s.PC) { MAXMSG("MeuAnalysis",Msg::kInfo,100) <<"Entry="<<entry<<", run="<<rec.GetRun() <<", snarl="<<rec.GetSnarl()<<endl; snarlList<<s.Run<<"\t"<<s.Snarl<<endl; } }//end of for MSG("MeuAnalysis",Msg::kInfo)<<"Finished main loop"<<endl; MSG("MeuAnalysis",Msg::kInfo) <<" ** Finished SnarlList method **"<<endl; }

void MeuAnalysis::SpillPlots (   )

<<PCCounter-bothSMCounter-coilHitCounter- Definition at line 3273 of file MeuAnalysis.cxx. References abs(), NtpSRPlane::beg, MeuReco::CalcFCPC(), MeuReco::CalcStripDists(), MeuReco::CalcVtx(), MeuReco::CalcVtxSpecialND(), MeuReco::CalcWindow(), MeuCuts::CheckTrkViews(), MeuReco::CheckWinContainment(), MeuSummary::Clear(), MeuSummary::DistToEdgeFid, NtpSRPlane::end, MeuSummary::EndDistToEdge, MeuSummary::EndPlane, MeuSummary::EndY, MeuSummary::Event, NtpStRecord::evt, MeuCuts::ExtractPlInfo(), MeuCuts::ExtractTruthInfo(), MeuSummary::FC, MeuCuts::FillEnVsDist(), MeuCuts::FillEnVsDist2(), MeuCuts::FillSTSumDetails(), MeuCuts::FillSTSumRecoDetails(), MeuCuts::FillTimeHistos(), MeuCuts::FilterBadDistEndStrip(), MeuCuts::FilterBadEndGaps(), MeuCuts::FilterBadEvtPerSlc(), MeuCuts::FilterBadTrackEnd(), MeuCuts::FilterBadTrkTimes(), fOutFile, fRun, RecRecordImp< T >::GetHeader(), RecDataHeader::GetRun(), RecPhysicsHeader::GetSnarl(), InitialiseLoopVariables(), MeuHitInfo::LPos, MAXMSG, MeuSummary::MCEndPlane, MeuSummary::MCHighEn, MeuSummary::MCLowEn, MeuSummary::MCVtxPlane, MSG, NtpSREvent::nshower, NtpSRTrack::nstrip, NtpSREvent::ntrack, OpenFile(), OpenTxtFile(), MeuSummary::PC, NtpSRShower::plane, MeuHitInfo::Plane, NtpSRTrack::plane, print(), MeuCuts::PrintNtpSt(), Recalibrate(), MeuReco::Reconstruct(), MeuSummary::RFid, MeuSummary::Run, s(), SetLoopVariables(), NtpSREvent::shw, NtpStRecord::shw, MeuHitInfo::SigMap, MeuSummary::SimFlag, NtpSREvent::slc, MeuSummary::Snarl, MeuHitInfo::TPos, MeuSummary::TrigSrc, NtpSREvent::trk, NtpStRecord::trk, MeuSummary::VtxDistToEdge, MeuSummary::VtxPlane, MeuSummary::VtxY, MeuSummary::WinAdc, MeuSummary::WinPe, MeuSummary::WinSigCor, MeuSummary::WinSigLin, MeuSummary::WinSigMap, MeuSummary::WinVtxSidePl, MeuHitInfo::X, MeuHitInfo::Y, and MeuHitInfo::Z. { MSG("MeuAnalysis",Msg::kInfo) <<" ** Running SpillPlots method... **"<<endl; //open the output file for the histograms fOutFile=this->OpenFile(100,"SpillPlots"); TH1F* hWinSigMap=new TH1F("hWinSigMap","hWinSigMap",3000,-1,4000); hWinSigMap->SetTitle("hWinSigMap"); hWinSigMap->GetXaxis()->SetTitle("hWinSigMap"); hWinSigMap->GetXaxis()->CenterTitle(); hWinSigMap->SetFillColor(0); hWinSigMap->SetLineWidth(3); //hWinSigMap->SetBit(TH1::kCanRebin); TH1F* hWinSigCor=new TH1F("hWinSigCor","hWinSigCor",3000,-1,4000); hWinSigCor->SetTitle("hWinSigCor"); hWinSigCor->GetXaxis()->SetTitle("hWinSigCor"); hWinSigCor->GetXaxis()->CenterTitle(); hWinSigCor->SetFillColor(0); hWinSigCor->SetLineWidth(3); //hWinSigCor->SetBit(TH1::kCanRebin); TH1F* hWinSigLin=new TH1F("hWinSigLin","hWinSigLin",3000,-1,4000); hWinSigLin->SetTitle("hWinSigLin"); hWinSigLin->GetXaxis()->SetTitle("hWinSigLin"); hWinSigLin->GetXaxis()->CenterTitle(); hWinSigLin->SetFillColor(0); hWinSigLin->SetLineWidth(3); //hWinSigLin->SetBit(TH1::kCanRebin); TH1F* hWinAdc=new TH1F("hWinAdc","hWinAdc",3000,-1,4000); hWinAdc->SetTitle("hWinAdc"); hWinAdc->GetXaxis()->SetTitle("hWinAdc"); hWinAdc->GetXaxis()->CenterTitle(); hWinAdc->SetFillColor(0); hWinAdc->SetLineWidth(3); //hWinAdc->SetBit(TH1::kCanRebin); TH1F* hWinPe=new TH1F("hWinPe","hWinPe",3000,-1,400); hWinPe->SetTitle("hWinPe"); hWinPe->GetXaxis()->SetTitle("hWinPe"); hWinPe->GetXaxis()->CenterTitle(); hWinPe->SetFillColor(0); hWinPe->SetLineWidth(3); //hWinPe->SetBit(TH1::kCanRebin); TH1F* hTrkVtxPl=new TH1F("hTrkVtxPl","hTrkVtxPl",501,-1,500); hTrkVtxPl->SetTitle("hTrkVtxPl"); hTrkVtxPl->GetXaxis()->SetTitle("hTrkVtxPl"); hTrkVtxPl->GetXaxis()->CenterTitle(); hTrkVtxPl->SetFillColor(0); hTrkVtxPl->SetLineWidth(3); //hTrkVtxPl->SetBit(TH1::kCanRebin); TH1F* hTrkEndPl=new TH1F("hTrkEndPl","hTrkEndPl",501,-1,500); hTrkEndPl->SetTitle("hTrkEndPl"); hTrkEndPl->GetXaxis()->SetTitle("hTrkEndPl"); hTrkEndPl->GetXaxis()->CenterTitle(); hTrkEndPl->SetFillColor(0); hTrkEndPl->SetLineWidth(3); //hTrkEndPl->SetBit(TH1::kCanRebin); TH1F* hTrkLengthPl=new TH1F("hTrkLengthPl","hTrkLengthPl",501,-1,500); hTrkLengthPl->SetTitle("hTrkLengthPl"); hTrkLengthPl->GetXaxis()->SetTitle("hTrkLengthPl"); hTrkLengthPl->GetXaxis()->CenterTitle(); hTrkLengthPl->SetFillColor(0); hTrkLengthPl->SetLineWidth(3); //hTrkLengthPl->SetBit(TH1::kCanRebin); TH1F* hWinPlFromShwEnd=new TH1F ("hWinPlFromShwEnd","hWinPlFromShwEnd",1000,-500,500); hWinPlFromShwEnd->SetTitle("hWinPlFromShwEnd"); hWinPlFromShwEnd->GetXaxis()->SetTitle("hWinPlFromShwEnd"); hWinPlFromShwEnd->GetXaxis()->CenterTitle(); hWinPlFromShwEnd->SetFillColor(0); hWinPlFromShwEnd->SetLineWidth(3); //hWinPlFromShwEnd->SetBit(TH1::kCanRebin); TH1F* hShwVtxPl=new TH1F("hShwVtxPl","hShwVtxPl",501,-1,500); hShwVtxPl->SetTitle("hShwVtxPl"); hShwVtxPl->GetXaxis()->SetTitle("hShwVtxPl"); hShwVtxPl->GetXaxis()->CenterTitle(); hShwVtxPl->SetFillColor(0); hShwVtxPl->SetLineWidth(3); //hShwVtxPl->SetBit(TH1::kCanRebin); TH1F* hShwEndPl=new TH1F("hShwEndPl","hShwEndPl",501,-1,500); hShwEndPl->SetTitle("hShwEndPl"); hShwEndPl->GetXaxis()->SetTitle("hShwEndPl"); hShwEndPl->GetXaxis()->CenterTitle(); hShwEndPl->SetFillColor(0); hShwEndPl->SetLineWidth(3); //hShwEndPl->SetBit(TH1::kCanRebin); TH1F* hShwLengthPl=new TH1F("hShwLengthPl","hShwLengthPl",501,-1,500); hShwLengthPl->SetTitle("hShwLengthPl"); hShwLengthPl->GetXaxis()->SetTitle("hShwLengthPl"); hShwLengthPl->GetXaxis()->CenterTitle(); hShwLengthPl->SetFillColor(0); hShwLengthPl->SetLineWidth(3); //hShwLengthPl->SetBit(TH1::kCanRebin); TH1F* hTrkShwVtxDiff=new TH1F("hTrkShwVtxDiff","hTrkShwVtxDiff", 1000,-500,500); hTrkShwVtxDiff->SetTitle("hTrkShwVtxDiff"); hTrkShwVtxDiff->GetXaxis()->SetTitle("hTrkShwVtxDiff"); hTrkShwVtxDiff->GetXaxis()->CenterTitle(); hTrkShwVtxDiff->SetFillColor(0); hTrkShwVtxDiff->SetLineWidth(3); //hTrkShwVtxDiff->SetBit(TH1::kCanRebin); TH1F* hTrkShwEndDiff=new TH1F("hTrkShwEndDiff","hTrkShwEndDiff", 1000,-500,500); hTrkShwEndDiff->SetTitle("hTrkShwEndDiff"); hTrkShwEndDiff->GetXaxis()->SetTitle("hTrkShwEndDiff"); hTrkShwEndDiff->GetXaxis()->CenterTitle(); hTrkShwEndDiff->SetFillColor(0); hTrkShwEndDiff->SetLineWidth(3); //hTrkShwEndDiff->SetBit(TH1::kCanRebin); TH1F* hShwVtxPl2=new TH1F("hShwVtxPl2","hShwVtxPl2",501,-1,500); hShwVtxPl2->SetTitle("hShwVtxPl2"); hShwVtxPl2->GetXaxis()->SetTitle("hShwVtxPl2"); hShwVtxPl2->GetXaxis()->CenterTitle(); hShwVtxPl2->SetFillColor(0); hShwVtxPl2->SetLineWidth(3); //hShwVtxPl2->SetBit(TH1::kCanRebin); TH1F* hShwEndPl2=new TH1F("hShwEndPl2","hShwEndPl2",501,-1,500); hShwEndPl2->SetTitle("hShwEndPl2"); hShwEndPl2->GetXaxis()->SetTitle("hShwEndPl2"); hShwEndPl2->GetXaxis()->CenterTitle(); hShwEndPl2->SetFillColor(0); hShwEndPl2->SetLineWidth(3); //hShwEndPl2->SetBit(TH1::kCanRebin); TH1F* hShwLengthPl2=new TH1F("hShwLengthPl2","hShwLengthPl2",501,-1,500); hShwLengthPl2->SetTitle("hShwLengthPl2"); hShwLengthPl2->GetXaxis()->SetTitle("hShwLengthPl2"); hShwLengthPl2->GetXaxis()->CenterTitle(); hShwLengthPl2->SetFillColor(0); hShwLengthPl2->SetLineWidth(3); //hShwLengthPl2->SetBit(TH1::kCanRebin); TH1F* hTrkShwVtxDiff2=new TH1F("hTrkShwVtxDiff2","hTrkShwVtxDiff2", 1000,-500,500); hTrkShwVtxDiff2->SetTitle("hTrkShwVtxDiff2"); hTrkShwVtxDiff2->GetXaxis()->SetTitle("hTrkShwVtxDiff2"); hTrkShwVtxDiff2->GetXaxis()->CenterTitle(); hTrkShwVtxDiff2->SetFillColor(0); hTrkShwVtxDiff2->SetLineWidth(3); //hTrkShwVtxDiff2->SetBit(TH1::kCanRebin); TH1F* hTrkShwEndDiff2=new TH1F("hTrkShwEndDiff2","hTrkShwEndDiff2", 1000,-500,500); hTrkShwEndDiff2->SetTitle("hTrkShwEndDiff2"); hTrkShwEndDiff2->GetXaxis()->SetTitle("hTrkShwEndDiff2"); hTrkShwEndDiff2->GetXaxis()->CenterTitle(); hTrkShwEndDiff2->SetFillColor(0); hTrkShwEndDiff2->SetLineWidth(3); //hTrkShwEndDiff2->SetBit(TH1::kCanRebin); TH1F* hDiffYEnds=new TH1F("hDiffYEnds","hDiffYEnds",201,-4.5,4.5); hDiffYEnds->SetTitle("hDiffYEnds"); hDiffYEnds->GetXaxis()->SetTitle("hDiffYEnds"); hDiffYEnds->GetXaxis()->CenterTitle(); hDiffYEnds->SetFillColor(0); hDiffYEnds->SetLineWidth(3); //hDiffYEnds->SetBit(TH1::kCanRebin); TH1F* hDiffYEndsSR=new TH1F("hDiffYEndsSR","hDiffYEndsSR",201,-4.5,4.5); hDiffYEndsSR->SetTitle("hDiffYEndsSR"); hDiffYEndsSR->GetXaxis()->SetTitle("hDiffYEndsSR"); hDiffYEndsSR->GetXaxis()->CenterTitle(); hDiffYEndsSR->SetFillColor(0); hDiffYEndsSR->SetLineWidth(3); //hDiffYEndsSR->SetBit(TH1::kCanRebin); TH1F* hEvtShws=new TH1F("hEvtShws","hEvtShws",100,-50,50); hEvtShws->SetTitle("hEvtShws"); hEvtShws->GetXaxis()->SetTitle("hEvtShws"); hEvtShws->GetXaxis()->CenterTitle(); hEvtShws->SetFillColor(0); hEvtShws->SetLineWidth(3); //hEvtShws->SetBit(TH1::kCanRebin); TProfile* pEnVsPl=new TProfile("pEnVsPl","pEnVsPl",1000,-500,500); pEnVsPl->SetTitle("Signal vs Distance from Trk Vtx"); pEnVsPl->GetXaxis()->SetTitle("Distance (Planes)"); pEnVsPl->GetXaxis()->CenterTitle(); pEnVsPl->GetYaxis()->SetTitle("Signal"); pEnVsPl->GetYaxis()->CenterTitle(); pEnVsPl->SetLineColor(1); pEnVsPl->SetFillColor(0); //pEnVsPl->SetBit(TH1::kCanRebin); TProfile* pEnVsPlShw=new TProfile ("pEnVsPlShw","pEnVsPlShw",1000,-500,500); pEnVsPlShw->SetTitle("Signal vs Distance from Trk Vtx (w/ shw)"); pEnVsPlShw->GetXaxis()->SetTitle("Distance (Planes)"); pEnVsPlShw->GetXaxis()->CenterTitle(); pEnVsPlShw->GetYaxis()->SetTitle("Signal"); pEnVsPlShw->GetYaxis()->CenterTitle(); pEnVsPlShw->SetLineColor(1); pEnVsPlShw->SetFillColor(0); //pEnVsPlShw->SetBit(TH1::kCanRebin); TProfile* pEnVsPlNoShw=new TProfile ("pEnVsPlNoShw","pEnVsPlNoShw",1000,-500,500); pEnVsPlNoShw->SetTitle("Signal vs Distance from Trk Vtx (no shw)"); pEnVsPlNoShw->GetXaxis()->SetTitle("Distance (Planes)"); pEnVsPlNoShw->GetXaxis()->CenterTitle(); pEnVsPlNoShw->GetYaxis()->SetTitle("Signal"); pEnVsPlNoShw->GetYaxis()->CenterTitle(); pEnVsPlNoShw->SetLineColor(1); pEnVsPlNoShw->SetFillColor(0); //pEnVsPlNoShw->SetBit(TH1::kCanRebin); TProfile* pEnVsPlFromShw=new TProfile ("pEnVsPlFromShw","pEnVsPlFromShw",1000,-500,500); pEnVsPlFromShw->SetTitle("Signal vs Distance from Shw End"); pEnVsPlFromShw->GetXaxis()->SetTitle("Distance (Planes)"); pEnVsPlFromShw->GetXaxis()->CenterTitle(); pEnVsPlFromShw->GetYaxis()->SetTitle("Signal"); pEnVsPlFromShw->GetYaxis()->CenterTitle(); pEnVsPlFromShw->SetLineColor(1); pEnVsPlFromShw->SetFillColor(0); //pEnVsPlFromShw->SetBit(TH1::kCanRebin); Int_t evtCounter=0; Int_t oneTrackCounter=0; Int_t oneEvtPerSlcCounter=0; Int_t longTrackCounter=0; Int_t closeTimesCounter=0; Int_t endPlaneSpectCounter=0; Int_t longTrack2Counter=0; Int_t badDirCounter=0; Int_t badDirSRCounter=0; Int_t FCCounter=0; Int_t PCCounter=0; Int_t TGCounter=0; Int_t badWindowCounter=0; Int_t badRecoCounter=0; Int_t goodWindowCounter=0; Int_t badFidCounter=0; Int_t badViewsCounter=0; Int_t badEndGapsCounter=0; Int_t badTrackEndCounter=0; Int_t bad2ndContCounter=0; Int_t badDistEndCounter=0; Int_t badShwDistCounter=0; string temps="myeventsLowEn"; if (fRun>10000) temps+="MC"; ofstream& eventsTxtFile=*(this->OpenTxtFile(100,"myevents")); ofstream& txtFileLowEn=*(this->OpenTxtFile(100,temps.c_str())); ofstream& txtFileNonPC=*(this->OpenTxtFile(100,"myeventsNonPC")); //reconstruction object MeuReco reco; MeuCuts cuts; MeuSummary s; //variable to turn on all the useful messages if required Msg::LogLevel_t logLevel=Msg::kDebug; this->InitialiseLoopVariables(); //for(Int_t entry=55000;entry<fEntries;entry++){ for(Int_t entry=0;entry<fEntries;entry++){ this->SetLoopVariables(entry); MSG("MeuAnalysis",logLevel) <<"Entry="<<entry<<endl; //make a reference instead of a pointer NtpStRecord& ntp=(*fRec); const RecCandHeader& rec=ntp.GetHeader(); MAXMSG("MeuAnalysis",Msg::kInfo,1) <<"First entry: run="<<rec.GetRun() <<", snarl="<<rec.GetSnarl()<<endl; TClonesArray& evtTca=(*ntp.evt); const Int_t numEvts=evtTca.GetEntriesFast(); //loop over the events for (Int_t ievt=0;ievt<numEvts;++ievt) { const NtpSREvent* pevt= dynamic_cast<NtpSREvent*>(evtTca[ievt]); const NtpSREvent& evt=(*pevt); evtCounter++; //ensure there is only 1 track in the event if (evt.ntrack!=1) continue; oneTrackCounter++; //filter out events occuring in a slice with more than 1 event if (cuts.FilterBadEvtPerSlc(ntp,evt.slc,entry)) continue; oneEvtPerSlcCounter++; TClonesArray& trkTca=(*ntp.trk); const NtpSRTrack* ptrk= dynamic_cast<NtpSRTrack*>(trkTca[evt.trk[0]]); const NtpSRTrack& trk=(*ptrk); //MeuSummary s; s.Clear(); s.RFid=1.0; s.DistToEdgeFid=0.4; s.Event=entry; cuts.FillSTSumDetails(ntp,s,ievt,evt.slc); //a map to hold the info for each plane in the entry map<Int_t,MeuHitInfo> plInfo; //cut on number of planes if (sqrt(pow(1.*trk.plane.beg-trk.plane.end,2))<8) continue; longTrackCounter++; MSG("MeuAnalysis",logLevel) <<"Entry="<<entry<<", run="<<rec.GetRun() <<", snarl="<<rec.GetSnarl() <<", trk.nstrips="<<trk.nstrip<<endl; if (s.Snarl==-1){ MSG("MeuAnalysis",Msg::kInfo) <<"Entry="<<entry<<", run="<<rec.GetRun() <<", snarl="<<rec.GetSnarl() <<", trk.nstrips="<<trk.nstrip<<endl; cuts.PrintNtpSt(ntp); } if (cuts.FilterBadTrkTimes(ntp,s,evt)) continue; closeTimesCounter++; cuts.ExtractPlInfo(ntp,s,plInfo,evt); reco.CalcVtx(s,plInfo); hDiffYEndsSR->Fill(plInfo[s.VtxPlane].Y- plInfo[s.EndPlane].Y); cuts.ExtractTruthInfo(ntp,s,plInfo); Bool_t mcPosDir=s.MCVtxPlane<s.MCEndPlane; Bool_t srPosDir=trk.plane.beg<trk.plane.end; Bool_t posDir=s.VtxPlane<s.EndPlane; if (posDir!=mcPosDir && s.SimFlag==SimFlag::kMC){ badDirCounter++; MAXMSG("MeuAnalysis",Msg::kVerbose,200) <<"Direction wrong! Trk: truth " <<s.MCVtxPlane<<" -> "<<s.MCEndPlane <<", reco "<<s.VtxPlane<<" -> "<<s.EndPlane //<<", y="<<s.VtxY<<" -> "<<s.EndY <<endl; } if (srPosDir!=mcPosDir && s.SimFlag==SimFlag::kMC){ badDirSRCounter++; MAXMSG("MeuAnalysis",Msg::kVerbose,100) <<"Direction wrong in SR! Trk: truth " <<s.MCVtxPlane<<" -> "<<s.MCEndPlane <<", SR reco "<<trk.plane.beg<<" -> "<<trk.plane.end <<endl; } //check if track stops in spectrometer if (s.EndPlane>120) continue; endPlaneSpectCounter++; MSG("MeuAnalysis",logLevel)<<"CalcFCPC..."<<endl; reco.CalcFCPC(s,plInfo); //recalculate the vtx of the tracks //coming in from the spectrometer MSG("MeuAnalysis",logLevel)<<"CalcVtxSpecialND..."<<endl; reco.CalcVtxSpecialND(s,plInfo); //cut on number of planes with potential new vertex if (sqrt(pow(1.*s.VtxPlane-s.EndPlane,2))<8) continue; longTrack2Counter++; Bool_t TG=!s.FC && !s.PC; if (s.PC) PCCounter++; else if (s.FC) FCCounter++; else if (TG) TGCounter++; if (s.FC && s.PC) MSG("MeuAnalysis",Msg::kError) <<"Both FC and PC!"<<endl; if (s.PC) { Int_t diffAtEnd=-1; Int_t diffAtVtx=-1; if (!cuts.CheckTrkViews(s,&diffAtVtx,&diffAtEnd)) { badViewsCounter++; continue; } MSG("MeuAnalysis",logLevel)<<"Found PC event..."<<endl; Bool_t goodReco=reco.Reconstruct(s,plInfo); if (goodReco){ MSG("MeuAnalysis",logLevel)<<"Recalibrating..."<<endl; this->Recalibrate(ntp,s,plInfo,evt); MSG("MeuAnalysis",logLevel)<<"CalcWindow..."<<endl; reco.CalcWindow(s,plInfo); //recalculate the containment! s.DistToEdgeFid=0.35;//allow a slight shift in position reco.CalcFCPC(s,plInfo); reco.CalcStripDists(s,plInfo); MSG("MeuAnalysis",logLevel)<<"CheckContainment..."<<endl; //fill the tree once event has been analysed if (s.WinSigMap<=0) { badWindowCounter++; } else if (!reco.CheckWinContainment(s,plInfo)){ badFidCounter++; //MSG("MeuAnalysis",Msg::kInfo)<<"Outside fid:"<<endl; } else if (cuts.FilterBadTrackEnd(s,plInfo)){ badTrackEndCounter++; } else if (cuts.FilterBadEndGaps(s,plInfo)){ badEndGapsCounter++; } else if (!s.PC){ bad2ndContCounter++; txtFileNonPC<<entry<<endl; Bool_t print=false; if (print){ MAXMSG("MeuAnalysis",Msg::kDebug,50) <<endl<<endl<<"Recalculated containment and not PC!" <<" PC="<<s.PC<<", FC="<<s.FC <<" entry="<<entry <<", snarl="<<s.Snarl<<", run="<<s.Run<<endl <<"Low En="<<s.MCLowEn<<" high="<<s.MCHighEn<<endl <<endl; MAXMSG("MeuAnalysis",Msg::kDebug,50) <<"vtx Pl="<<s.VtxPlane<<", end Pl="<<s.EndPlane <<", dVtx="<<s.VtxDistToEdge <<", dEnd="<<s.EndDistToEdge<<endl; for (map<Int_t,MeuHitInfo>::const_iterator it= plInfo.begin(); it!=plInfo.end();++it){ const MeuHitInfo& c=it->second; MAXMSG("MeuAnalysis",Msg::kInfo,500) <<"pl="<<c.Plane<<", z="<<c.Z <<", x="<<c.X<<", y="<<c.Y <<" t="<<c.TPos<<", l="<<c.LPos<<endl; } } } else if (cuts.FilterBadDistEndStrip(s,plInfo)){ badDistEndCounter++; } else { hTrkVtxPl->Fill(s.VtxPlane); hTrkEndPl->Fill(s.EndPlane); hTrkLengthPl->Fill(s.EndPlane-s.VtxPlane); hDiffYEnds->Fill(s.VtxY-s.EndY); hEvtShws->Fill(evt.nshower); Int_t shwVtxPl=999; Int_t shwEndPl=-999; //loop over showers in evt TClonesArray& shwTca=(*ntp.shw); for(int i=0;i<evt.nshower;i++) { const NtpSRShower* pshw= dynamic_cast<NtpSRShower*>(shwTca[evt.shw[i]]); const NtpSRShower& shw=(*pshw); hShwVtxPl->Fill(shw.plane.beg); hShwEndPl->Fill(shw.plane.end); hShwLengthPl->Fill(shw.plane.end-shw.plane.beg); hTrkShwVtxDiff->Fill(s.VtxPlane-shw.plane.beg); hTrkShwEndDiff->Fill(s.EndPlane-shw.plane.end); //define a vtx shower to be one that starts within 5 pl Bool_t vtxShw=abs(s.VtxPlane-shw.plane.beg)<=5; if (vtxShw){ if (shw.plane.beg<shwVtxPl) shwVtxPl=shw.plane.beg; if (shw.plane.end>shwEndPl) shwEndPl=shw.plane.end; } } Bool_t vtxShw=false; if (shwEndPl!=-999 && shwVtxPl!=999) { vtxShw=true; hShwVtxPl2->Fill(shwVtxPl); hShwEndPl2->Fill(shwEndPl); hShwLengthPl2->Fill(shwEndPl-shwVtxPl); hTrkShwVtxDiff2->Fill(s.VtxPlane-shwVtxPl); hTrkShwEndDiff2->Fill(s.EndPlane-shwEndPl); } //fill profiles typedef map<Int_t,MeuHitInfo>::const_iterator plIt; for (plIt it=plInfo.begin();it!=plInfo.end();++it){ const MeuHitInfo& cp=it->second; //check that plane is within track (assume forward going) if (cp.Plane<s.VtxPlane || cp.Plane>s.EndPlane) continue; Int_t plFromVtx=cp.Plane-s.VtxPlane; pEnVsPl->Fill(plFromVtx,cp.SigMap); if (vtxShw) { pEnVsPlShw->Fill(plFromVtx,cp.SigMap); Int_t plFromShwEnd=cp.Plane-shwEndPl; if (plFromShwEnd>=0){ pEnVsPlFromShw->Fill(plFromShwEnd,cp.SigMap); } } else pEnVsPlNoShw->Fill(plFromVtx,cp.SigMap); } if (vtxShw){ Int_t winPlFromShwEnd=s.WinVtxSidePl-shwEndPl; hWinPlFromShwEnd->Fill(winPlFromShwEnd); //require that shw stops 10 planes before window starts if (winPlFromShwEnd<10 && s.TrigSrc>100) { badShwDistCounter++; continue; } } goodWindowCounter++; cuts.ExtractTruthInfo(ntp,s,plInfo); cuts.FillSTSumRecoDetails(s,plInfo); hWinSigMap->Fill(s.WinSigMap); hWinSigCor->Fill(s.WinSigCor); hWinSigLin->Fill(s.WinSigLin); hWinAdc->Fill(s.WinAdc); hWinPe->Fill(s.WinPe); eventsTxtFile<<entry<<endl; MAXMSG("MeuAnalysis",Msg::kInfo,10) <<"Good window: run=" <<rec.GetRun()<<", snarl="<<rec.GetSnarl() <<endl; if (s.MCLowEn>0.5*Munits::GeV && s.SimFlag==SimFlag::kMC) { txtFileLowEn<<entry<<endl; MSG("MeuAnalysis",Msg::kInfo) <<"Low En="<<s.MCLowEn<<" high="<<s.MCHighEn <<", vtxPl="<<s.VtxPlane<<", endPl="<<s.EndPlane <<endl<<"entry="<<entry <<", snarl="<<s.Snarl<<", run="<<s.Run<<endl; } cuts.FillEnVsDist(s,plInfo); cuts.FillEnVsDist2(s,plInfo); cuts.FillTimeHistos(ntp,evt,s); } } else { badRecoCounter++; MSG("MeuAnalysis",logLevel)<<"Bad reco..."<<endl; } } }//end of loop over events }//end of for MSG("MeuAnalysis",Msg::kInfo)<<"Finished main loop"<<endl; Double_t quantile=0.5;//quantile to compute Double_t meuWinSigMap=-1; hWinSigMap->GetQuantiles(1,&meuWinSigMap,&quantile); Double_t errWinSigMap=hWinSigMap->GetRMS()/ sqrt(hWinSigMap->GetEntries()); Double_t relErrWinSigMap=errWinSigMap/meuWinSigMap; Double_t meuWinSigCor=-1; hWinSigCor->GetQuantiles(1,&meuWinSigCor,&quantile); Double_t errWinSigCor=hWinSigCor->GetRMS()/ sqrt(hWinSigCor->GetEntries()); Double_t relErrWinSigCor=errWinSigCor/meuWinSigCor; Double_t meuWinSigLin=-1; hWinSigLin->GetQuantiles(1,&meuWinSigLin,&quantile); Double_t errWinSigLin=hWinSigLin->GetRMS()/ sqrt(hWinSigLin->GetEntries()); Double_t relErrWinSigLin=errWinSigLin/meuWinSigLin; Double_t meuWinAdc=-1; hWinAdc->GetQuantiles(1,&meuWinAdc,&quantile); Double_t errWinAdc=hWinAdc->GetRMS()/ sqrt(hWinAdc->GetEntries()); Double_t relErrWinAdc=errWinAdc/meuWinAdc; Double_t meuWinPe=-1; hWinPe->GetQuantiles(1,&meuWinPe,&quantile); Double_t errWinPe=hWinPe->GetRMS()/ sqrt(hWinPe->GetEntries()); Double_t relErrWinPe=errWinPe/meuWinPe; hWinPe->Fit("gaus"); TF1* fPe=hWinPe->GetFunction("gaus"); Double_t errWinPe2=fPe->GetParameter(2)/sqrt(hWinPe->GetEntries()); // cout<<"RESOLUTIONS (from Gaussian fit):"<<endl // <<"ADCs: "<<fAdc->GetParameter(2)/fAdc->GetParameter(1)<<endl // <<"PEs: "<<fPe->GetParameter(2)/fPe->GetParameter(1)<<endl MSG("MeuAnalysis",Msg::kInfo) <<endl<<"*** Run Summary ***"<<endl //<<"Raw MEU="<<meu<<" +/- "<<err<<" ("<<100.*relErr<<"%)"<<endl <<"Raw MEU (WinSigMap)="<<meuWinSigMap<<" +/- "<<errWinSigMap <<" ("<<100.*relErrWinSigMap<<"%)"<<endl <<"Raw MEU (WinSigCor)="<<meuWinSigCor<<" +/- "<<errWinSigCor <<" ("<<100.*relErrWinSigCor<<"%)"<<endl <<"Raw MEU (WinSigLin)="<<meuWinSigLin<<" +/- "<<errWinSigLin <<" ("<<100.*relErrWinSigLin<<"%)"<<endl <<"Raw MEU (WinAdc) ="<<meuWinAdc<<" +/- "<<errWinAdc <<" ("<<100.*relErrWinAdc<<"%)"<<endl <<"Raw MEU (WinPe) ="<<meuWinPe<<" +/- "<<errWinPe <<" ("<<100.*relErrWinPe<<"%), fitErr="<<errWinPe2<<endl; //<<endl /*<<"Not LI="<<notLICounter<<endl <<"Not MM="<<notMMCounter<<endl <<"Good track="<<goodTrackCounter<<endl <<"Good length="<<goodLengthCounter<<endl <<"Good end & vtx="<<goodEndAndVtxCounter<<endl <<"Good XY reco="<<goodXYCounter<<endl <<endl <<"Good containmentPC="<<containmentCounter<<endl <<" both SM hit (when PC)="<<bothSMCounter <<" ("<<bothSM<<"% of PC)"<<endl <<" coil hit (when PC) ="<<coilHitCounter <<" ("<<coil_PC<<"% of PC)"<<endl <<" close LI (when PC) ="<<closeLICounter <<" ("<<closeLI<<"% of PC)"<<endl <<" dead Chips (when PC) ="<<deadChipsCounter <<" ("<<deadChips<<"% of PC)"<<endl <<" boundary TF (when PC)="<<boundaryTFCounter <<" ("<<boundaryTF<<"% of PC)"<<endl <<" bad views (when PC) ="<<badViewsCounter <<" ("<<badViews<<"% of PC)"<<endl <<" bad trk end (when PC)="<<badTrkEndCounter <<" ("<<badTrkEnd<<"% of PC)"<<endl <<" bad trk gap (when PC)="<<badTrkGapsCounter <<" ("<<badTrkGaps<<"% of PC)"<<endl <<" bad traceZ (when PC) ="<<badTraceZCounter <<" ("<<badTraceZ<<"% of PC)"<<endl <<" bad window (when PC) ="<<badWindowCounter <<" ("<<badWin<<"% of PC)"<<endl <<"Good containmentFC="<<containmentCounterFC<<endl <<" both SM hit (when FC)="<<bothSMCounterFC <<" ("<<bothSMFC<<"% of FC)"<<endl <<" coil hit (when FC) ="<<coilHitCounterFC <<" ("<<coil_FC<<"% of FC)"<<endl <<"Good containmentTG="<<containmentCounterTG<<endl <<" both SM hit (when TG)="<<bothSMCounterTG <<" ("<<bothSMTG<<"% of TG)"<<endl <<" coil hit (when TG) ="<<coilHitCounterTG <<" ("<<coil_TG<<"% of TG)"<<endl <<endl <<"Bad charge pos moment="<<momentCounter<<endl <<"Too high StripsPerPlane="<<stPerPlCounter<<endl <<" Either of above two="<<stPlOrMomentCounter<<endl <<endl <<"MC Info:"<<endl <<" Muon lowest energy >0.6 GeV="<<lowEnCounter<<endl <<" Muon lowest energy >5.0 GeV="<<lowEn5Counter<<endl <<endl */ Int_t totalBadPCEvents=badWindowCounter+badFidCounter+ badRecoCounter+badViewsCounter+ badTrackEndCounter+badEndGapsCounter+bad2ndContCounter+ badDistEndCounter+badShwDistCounter; MSG("MeuAnalysis",Msg::kInfo) <<"Total Entries (snarls)="<<fEntries<<endl <<"Total Evts="<<evtCounter<<endl <<" One track ="<<oneTrackCounter<<endl <<" One evt in slc ="<<oneEvtPerSlcCounter<<endl <<" Long track ="<<longTrackCounter<<endl <<" Close trk times="<<closeTimesCounter<<endl <<" End plane cal ="<<endPlaneSpectCounter<<endl <<" Long track 2 ="<<longTrack2Counter<<endl <<" Bad dir (MC) ="<<badDirCounter<<endl <<" Bad dir SR (MC)="<<badDirSRCounter<<endl <<"Containment summary:"<<endl <<" PC="<<PCCounter<<endl // <<" (PC - all cuts=" //closeLICounter-deadChipsCounter-boundaryTFCounter- //badViewsCounter-badTrkEndCounter-badTrkGapsCounter-badTraceZCounter- //badWindowCounter<<")"<<endl <<" FC="<<FCCounter<<endl//<<" (FC-coil-bothSM=" //<<FCCounter-coilHitCounterFC-bothSMCounterFC<<")"<<endl <<" TG="<<TGCounter<<endl//<<" (TG-coil-bothSM=" // <<TGCounter-coilHitCounterTG-bothSMCounterTG<<")"<<endl <<" FC+PC+TG="<<FCCounter+PCCounter+TGCounter<<endl <<endl <<"Total PC with good track window = "<<goodWindowCounter<<endl <<" bad track window = "<<badWindowCounter<<endl <<" bad fiducial vol = "<<badFidCounter<<endl <<" bad reco = "<<badRecoCounter<<endl <<" bad views = "<<badViewsCounter<<endl <<" bad track end = "<<badTrackEndCounter<<endl <<" bad track end gaps= "<<badEndGapsCounter<<endl <<" bad 2nd contain. = "<<bad2ndContCounter<<endl <<" bad dist strip end= "<<badDistEndCounter<<endl <<" bad dist to shw = "<<badShwDistCounter<<endl <<" Total bad = "<<totalBadPCEvents<<endl <<" Total bad + good = " <<totalBadPCEvents+goodWindowCounter<<endl <<endl; if (fRec) delete fRec; MSG("MeuAnalysis",Msg::kInfo) <<" ** Finished SpillPlots method **"<<endl; }

void MeuAnalysis::StoreOrFinishSummaryTree (  MeuSummaryWriter *  psOut, MeuCounter *  pcnt, Bool_t  finish )  const

Definition at line 5018 of file MeuAnalysis.cxx. References MSG, MeuHistos::PrintMeuValues(), MeuCounter::PrintNtpSt(), and MeuSummaryWriter::SummaryTreeFinish(). Referenced by MeuCalModule::EndJob(), and MakeSummaryTreeWithNtpSt(). { //this stores the pointers the first time the function is called static MeuCounter& cnt=*pcnt; static MeuSummaryWriter& summaryOut=*psOut; //this section runs methods on objects at the end of job if (finish){ MSG("MeuAnalysis",Msg::kInfo) <<endl<<"*** Run Summary ***"<<endl; MeuHistos histos; histos.PrintMeuValues();//this looks up TObjects by name //print out the counter information cnt.PrintNtpSt(); //finish the summary tree and output it summaryOut.SummaryTreeFinish(); } }

void MeuAnalysis::Test (   )

Int_t ns = trk.nstrip - nsloop - 1; Definition at line 771 of file MeuAnalysis.cxx. References NtpSRFiducial::dr, NtpSRTrack::ds, NtpSRFiducial::dz, NtpSRTrack::end, NtpStRecord::evthdr, fChain, NtpSRTrack::fidend, fRec, MAXMSG, MSG, NtpSREventSummary::ndigit, NtpSRTrack::nstrip, NtpSREventSummary::ntrack, NtpSRStrip::ph1, NtpSRStrip::plane, NtpSRTrack::range, NtpSRPulseHeight::sigcor, NtpSRTrack::stp, NtpStRecord::stp, NtpSRTrack::stpz, NtpSRStrip::strip, NtpSRStrip::tpos, NtpStRecord::trk, NtpSRTrack::vtx, and NtpSRVertex::z. Referenced by OpenFile(). { MSG("MeuAnalysis",Msg::kInfo) <<" ** Running Test method... **"<<endl; //this->InitialiseLoopVariables(); // for(Int_t event=0;event<fEntries;event++){ //this->SetLoopVariables(event); for (int i=0;i<5;i++) { fChain->GetEntry(i); NtpStRecord& ntpst=(*fRec); // Make a reference instead of a pointer TClonesArray& tcaTk=(*ntpst.trk); Int_t numTrks=tcaTk.GetEntries(); MSG("MeuAnalysis",Msg::kInfo) <<"numTrks="<<numTrks <<", ndigit="<<fRec->evthdr.ndigit <<", nTrack="<<fRec->evthdr.ntrack<<endl; //Loop over tracks for (Int_t itrk=0;itrk<numTrks;itrk++){ const NtpSRTrack* ptrk= dynamic_cast<NtpSRTrack*>(tcaTk[itrk]); const NtpSRTrack& trk=(*ptrk); cout<<"range="<<trk.range<<" nstrip="<<trk.nstrip<<" ds="<<trk.ds <<" fidend.dr,dz="<<trk.fidend.dr<<" "<<trk.fidend.dz<<" end-vtx.z:" <<trk.end.z-trk.vtx.z<<" "<<trk.stpz[1]-trk.stpz[0]<<endl; TClonesArray& tcaStp=(*ntpst.stp); Int_t numStps=tcaStp.GetEntries(); MSG("MeuAnalysis",Msg::kInfo) <<"numStps="<<numStps<<endl; for (Int_t i=0;i<trk.nstrip;i++){ //check for bug where strip index is -1 if (trk.stp[i]<0) { MAXMSG("MeuAnalysis",Msg::kInfo,500) <<"Skipping strip with trk.stp[i]="<<trk.stp[i]<<endl; continue; } //const NtpSRStrip* pstp=dynamic_cast<NtpSRStrip*>(tcaStp[i]); const NtpSRStrip* pstp=dynamic_cast<NtpSRStrip*>(tcaStp[trk.stp[i]]); const NtpSRStrip& stp=(*pstp); MSG("MeuAnalysis",Msg::kInfo) <<"Strip="<<stp.strip<<", pl="<<stp.plane <<", sigCor="<<stp.ph1.sigcor<<", tpos="<<stp.tpos<<endl; } } TClonesArray& tcaStp=(*ntpst.stp); Int_t numStps=tcaStp.GetEntries(); MSG("MeuAnalysis",Msg::kInfo) <<endl<<"2nd: numStps="<<numStps<<endl; for (Int_t i=0;i<numStps;i++){ const NtpSRStrip* pstp=dynamic_cast<NtpSRStrip*>(tcaStp[i]); const NtpSRStrip& stp=(*pstp); MSG("MeuAnalysis",Msg::kInfo) <<"Strip="<<stp.strip<<", pl="<<stp.plane<<endl; } // for (Int_t nsloop=0;nsloop<trk.nstrip;nsloop++) { //const NtpSRStrip* pstp= //dynamic_cast<NtpSRStrip*>(tcaStp[trk.stp[ns]]); //const NtpSRStrip& stp=(*pstp); // if (stp.plane) ; // Avoid compiler message ifndef PRINTLOTS //float du = trk.stpu[ns]-uprev; //float dv = trk.stpv[ns]-vprev; //float dz = trk.stpz[ns]-zprev; //cout<<" "<<ns<<" "<<du<<" "<<dv<<" "<<dz // <<" uv:"<<trk.stpu[ns] <<" "<<trk.stpv[ns] // <<" xyz:"<<trk.stpx[ns] <<" "<<trk.stpy[ns]<<" "<<trk.stpz[ns] // <<" ds:"<<trk.stpds[ns]<<" stp: "<<trk.stp[ns] // <<" "<<stp.plane<<endl; // printf("%2d%8.4f%8.4f%8.4f%8.4f%2d uv:%8.4f%8.4f" // " xyz:%8.4f%8.4f%8.4f%8.4f" // " stp: %4d%4d\n" // ,ns,du,dv,dz,stp.tpos,stp.planeview // ,trk.stpu[ns],trk.stpv[ns] // ,trk.stpx[ns],trk.stpy[ns],trk.stpz[ns] // ,trk.stpds[ns],trk.stp[ns],stp.plane); //uprev = trk.stpu[ns]; //vprev = trk.stpv[ns]; //zprev = trk.stpz[ns]; //}// stp loop }//end of for MSG("MeuAnalysis",Msg::kInfo)<<"Finished main loop"<<endl; MSG("MeuAnalysis",Msg::kInfo) <<" ** Finished Test method **"<<endl; }

void MeuAnalysis::TestEventLoop (   )

Definition at line 888 of file MeuAnalysis.cxx. References NtpSREvent::end, NtpStRecord::evt, MeuCuts::FilterBadEvtPerSlc(), InitialiseLoopVariables(), MAXMSG, MSG, NtpSREvent::nshower, NtpSREvent::nstrip, NtpSRShower::nstrip, NtpSRTrack::nstrip, NtpSRSlice::nstrip, NtpSREvent::ntrack, NtpSRStrip::ph0, NtpSRStrip::ph1, NtpSRPulseHeight::raw, SetLoopVariables(), NtpSREvent::shw, NtpStRecord::shw, NtpSREvent::slc, NtpStRecord::slc, NtpSREvent::stp, NtpSRTrack::stp, NtpSRSlice::stp, NtpStRecord::stp, NtpSRStrip::time0, NtpSREvent::trk, and NtpStRecord::trk. { MSG("MeuAnalysis",Msg::kInfo) <<" ** Running TestEventLoop method... **"<<endl; MeuCuts cuts; this->InitialiseLoopVariables(); for(Int_t entry=0;entry<fEntries;entry++) { //for(Int_t entry=0;entry<5000;entry++) { this->SetLoopVariables(entry); NtpStRecord& ntp=(*fRec); TClonesArray& evtTca=(*ntp.evt); const Int_t numEvts=evtTca.GetEntriesFast(); TClonesArray& trkTca=(*ntp.trk); const Int_t numTrks=trkTca.GetEntriesFast(); TClonesArray& shwTca=(*ntp.shw); const Int_t numShws=shwTca.GetEntriesFast(); TClonesArray& stpTca=(*ntp.stp); const Int_t numStps=stpTca.GetEntries(); TClonesArray& slcTca=(*ntp.slc); const Int_t numSlcs=slcTca.GetEntries(); Float_t totalSnarlAdc=0; Float_t totalEvtAdc=0; Float_t totalSlcAdc=0; Float_t thisSlcAdc=0; for (Int_t i=0;i<numStps;i++) { const NtpSRStrip* pstp= dynamic_cast<NtpSRStrip*>(stpTca[i]); const NtpSRStrip& stp=(*pstp); totalSnarlAdc+=stp.ph0.raw+stp.ph1.raw; } if (numEvts==0 && numTrks>0) { cout<<"Ahhh, num events in tca="<<numEvts <<", trks="<<numTrks<<endl; } if (numEvts>1) { MAXMSG("MeuAnalysis",Msg::kInfo,10) <<"High num evts="<<numEvts<<", trks="<<numTrks<<endl; } if (numSlcs!=1) { MAXMSG("MeuAnalysis",Msg::kInfo,10) <<"High num slices in tca="<<numSlcs <<", trks="<<numTrks<<endl; //continue; } MAXMSG("MeuAnalysis",Msg::kInfo,100) <<"Num events in tca="<<numEvts <<", slcs="<<numSlcs<<", trks="<<numTrks<<endl; Int_t slcStps=0; //loop over the slices for (Int_t islc=0;islc<numSlcs;++islc) { const NtpSRSlice* pslc= dynamic_cast<NtpSRSlice*>(slcTca[islc]); const NtpSRSlice& slc=(*pslc); MAXMSG("MeuAnalysis",Msg::kDebug,200) <<" slice "<<islc+1<<" of "<<numSlcs<<endl; slcStps+=slc.nstrip; //loop over strips in slc for (Int_t i=0;i<slc.nstrip;++i) { const NtpSRStrip* pstp= dynamic_cast<NtpSRStrip*>(stpTca[slc.stp[i]]); const NtpSRStrip& stp=(*pstp); totalSlcAdc+=stp.ph0.raw+stp.ph1.raw; } //loop over tracks in slc ---- can't do this!!! //for(int i=0;i<slc.ntrack;i++){ //const NtpSRTrack* ptrk= // dynamic_cast<NtpSRTrack*>(trkTca[slc.trk[i]]); //const NtpSRTrack& trk=(*ptrk); //trkStps=trk.nstrip; //MAXMSG("MeuAnalysis",Msg::kInfo,300) // <<" track "<<i+1<<" of "<<slc.ntrack // <<", trkStps="<<trkStps<<endl; //} } map<Int_t,Int_t> evtPerSlc; //loop over the events for (Int_t ntpevt=0;ntpevt<numEvts;++ntpevt) { const NtpSREvent* pevt= dynamic_cast<NtpSREvent*>(evtTca[ntpevt]); const NtpSREvent& evt=(*pevt); MAXMSG("MeuAnalysis",Msg::kInfo,100) <<"Entry "<<entry<<" has tracks="<<evt.ntrack <<", shws="<<numShws<<", slcs="<<numSlcs<<endl; //get the slice associated with this event const NtpSRSlice* pslc= dynamic_cast<NtpSRSlice*>(slcTca[evt.slc]); const NtpSRSlice& slc=(*pslc); evtPerSlc[evt.slc]++; //check that the filter works Int_t evts=-1; Bool_t goodEvt=!cuts.FilterBadEvtPerSlc(ntp,evt.slc,entry,&evts); if (goodEvt){ MAXMSG("MeuAnalysis",Msg::kInfo,100) <<"TEST::entry="<<entry <<", evt="<<ntpevt<<" is good, evtsPerSlc="<<evts <<" (slc="<<evt.slc<<")"<<endl; } else{ MAXMSG("MeuAnalysis",Msg::kInfo,100) <<"TEST::entry="<<entry <<", evt="<<ntpevt<<" is bad, evtsPerSlc="<<evts <<" (slc="<<evt.slc<<")"<<endl; } if (evt.ntrack!=1 || evt.nshower>1) continue; Int_t trkStps=0; Int_t shwStps=0; //loop over strips in slc for (Int_t i=0;i<slc.nstrip;++i) { const NtpSRStrip* pstp= dynamic_cast<NtpSRStrip*>(stpTca[slc.stp[i]]); const NtpSRStrip& stp=(*pstp); thisSlcAdc+=stp.ph0.raw+stp.ph1.raw; } //loop over tracks in evt for(int i=0;i<evt.ntrack;i++) { const NtpSRTrack* ptrk= dynamic_cast<NtpSRTrack*>(trkTca[evt.trk[i]]); const NtpSRTrack& trk=(*ptrk); trkStps=trk.nstrip; for (Int_t i=0;i<trk.nstrip;i++) { //check for bug where strip index is -1 if (trk.stp[i]<0) { MAXMSG("MeuAnalysis",Msg::kInfo,500) <<"Skipping strip with trk.stp[i]="<<trk.stp[i]<<endl; continue; } const NtpSRStrip* pstp= dynamic_cast<NtpSRStrip*>(stpTca[trk.stp[i]]); const NtpSRStrip& stp=(*pstp); MAXMSG("MeuAnalysis",Msg::kDebug,100) <<"strip time="<<stp.time0<<endl; } } //loop over showers in evt for(int i=0;i<evt.nshower;i++) { const NtpSRShower* pshw= dynamic_cast<NtpSRShower*>(shwTca[evt.shw[i]]); const NtpSRShower& shw=(*pshw); shwStps=shw.nstrip; } //loop over strips in evt for (Int_t i=0;i<evt.nstrip;i++) { const NtpSRStrip* pstp= dynamic_cast<NtpSRStrip*>(stpTca[evt.stp[i]]); const NtpSRStrip& stp=(*pstp); totalEvtAdc+=stp.ph0.raw+stp.ph1.raw; } MAXMSG("MeuAnalysis",Msg::kInfo,100) <<"e="<<entry<<", evt="<<ntpevt <<", ADCs: snl="<<totalSnarlAdc <<", evt="<<totalEvtAdc <<", slc="<<totalSlcAdc <<", tslc="<<thisSlcAdc <<", #st="<<numStps<<","<<evt.nstrip<<","<<slcStps //<<", trk="<<trkStps<<", shw="<<shwStps <<endl; } //make this into a function //Bool_t FilterMultiEvtPerSlc(ntp,slc,entry) typedef map<Int_t,Int_t>::const_iterator slcIt; MAXMSG("MeuAnalysis",Msg::kInfo,30) <<"MAIN::Events per slice for entry="<<entry<<endl; for (slcIt it=evtPerSlc.begin();it!=evtPerSlc.end();++it) { MAXMSG("MeuAnalysis",Msg::kInfo,200) <<" slice "<<it->first<<" has "<<it->second<<" event(s)"<<endl; } }//end of for MSG("MeuAnalysis",Msg::kInfo)<<"Finished main loop"<<endl; MSG("MeuAnalysis",Msg::kInfo) <<" ** Finished TestEventLoop method **"<<endl; }

void MeuAnalysis::UseAtNu (   )  [static]

Definition at line 215 of file MeuAnalysis.cxx. { //fUseAtNu=true; }

void MeuAnalysis::WriteOutHistos (   )  const

Definition at line 222 of file MeuAnalysis.cxx. References fOutFile, and MSG. Referenced by MeuCalModule::EndJob(). { //write out the histos to the file, if it's open if (fOutFile){ if (fOutFile->IsWritable()){ fOutFile->cd(); //create a tree and set branches //TTree runInfo("runInfo","runInfo"); //have to create local variables here //I think it's something to do with the fact that they are //read in on other trees in this class as well??? //No, actually not because fSigCorsPerMip is not in a tree //don't understand this crazyness! //Int_t lSimFlag=0; //runInfo.Branch("SimFlag",&lSimFlag,"SimFlag/I"); //fill and write //runInfo.Fill(); //runInfo.Write(); MSG("MeuAnalysis",Msg::kInfo) <<"Writing histos to: "<<fOutFile->GetName()<<" ..."<<endl; fOutFile->Write(); //fOutFile->Close();//this makes histos disappear from canvases //so do it in the destructor (need to make MeuAnalysis on heap) } else { MSG("MeuAnalysis",Msg::kWarning) <<"File not writable!"<<endl; } } }

---

Member Data Documentation

TChain* MeuAnalysis::fChain [private]

Definition at line 101 of file MeuAnalysis.h. Referenced by MakeChain(), SetLoopVariables(), SetRootFileObjects(), and Test().

TChain* MeuAnalysis::fChainBD [private]

Definition at line 102 of file MeuAnalysis.h. Referenced by MakeChain(), SetLoopVariables(), and SetRootFileObjects().

Float_t MeuAnalysis::fEntries [private]

Definition at line 109 of file MeuAnalysis.h. Referenced by SetLoopVariables(), and SetRootFileObjects().

Float_t MeuAnalysis::fEntriesBD [private]

Definition at line 110 of file MeuAnalysis.h. Referenced by SetRootFileObjects().

AtmosEvent* MeuAnalysis::fev [private]

Definition at line 103 of file MeuAnalysis.h. Referenced by SetRootFileObjects().

string MeuAnalysis::fInputFileName = "" [static, private]

Definition at line 59 of file MeuAnalysis.cxx. Referenced by InputFileName(), and MakeFileList().

Bool_t MeuAnalysis::fLoadTrees = false [static, private]

Definition at line 60 of file MeuAnalysis.cxx. Referenced by LoadTrees().

TFile* MeuAnalysis::fOutFile [private]

Definition at line 104 of file MeuAnalysis.h. Referenced by BasicPlots(), N_1Plots(), SpillPlots(), WriteOutHistos(), and ~MeuAnalysis().

NtpStRecord* MeuAnalysis::fRec [private]

Definition at line 105 of file MeuAnalysis.h. Referenced by BasicPlots(), MakeSummaryTreeWithNtpSt(), SetRootFileObjects(), and Test().

Bool_t MeuAnalysis::fRecalibratePE = false [static, private]

Definition at line 63 of file MeuAnalysis.cxx. Referenced by RecalibratePE().

Bool_t MeuAnalysis::fRecalibrateSigCor = false [static, private]

Definition at line 62 of file MeuAnalysis.cxx. Referenced by RecalibrateSigCor().

Bool_t MeuAnalysis::fRecalibrateSigLin = false [static, private]

Definition at line 61 of file MeuAnalysis.cxx. Referenced by RecalibrateSigLin().

NtpBDLiteRecord* MeuAnalysis::fRecBD [private]

Definition at line 106 of file MeuAnalysis.h. Referenced by MakeSummaryTreeWithNtpSt(), and SetRootFileObjects().

Int_t MeuAnalysis::fRun [private]

Definition at line 111 of file MeuAnalysis.h. Referenced by MakeSummaryTreeWithAtNu(), and SpillPlots().

Int_t MeuAnalysis::fSubrun [private]

Definition at line 112 of file MeuAnalysis.h. Referenced by MakeSummaryTreeWithAtNu().

Bool_t MeuAnalysis::fUseAtNu [private]

Definition at line 107 of file MeuAnalysis.h. Referenced by MakeChain().

---

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

• MeuAnalysis.h
• MeuAnalysis.cxx

---