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

static void InputFileName (std::string f)
static void LoadTrees (Bool_t loadTrees)
static void UseAtNu ()
static void RecalibrateSigLin (Bool_t recalibrate=true)
static void RecalibrateSigCor (Bool_t recalibrate=true)
static 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
AtmosEventfev
TFile * fOutFile
NtpStRecordfRec
NtpBDLiteRecordfRecBD
Bool_t fUseAtNu
Float_t fEntries
Float_t fEntriesBD
Int_t fRun
Int_t fSubrun

Static Private Attributes

static std::string fInputFileName = ""
static Bool_t fLoadTrees = false
static Bool_t fRecalibrateSigLin = false
static Bool_t fRecalibrateSigCor = false
static Bool_t fRecalibratePE = false

Detailed Description

Definition at line 31 of file MeuAnalysis.h.


Constructor & Destructor Documentation

MeuAnalysis::MeuAnalysis (  ) 

Definition at line 71 of file MeuAnalysis.cxx.

References fChain, RecRecordImp< T >::GetHeader(), RecDataHeader::GetRun(), RecDataHeader::GetSubRun(), Msg::kInfo, and MSG.

00072 {
00073   MSG("MeuAnalysis",Msg::kInfo)
00074     <<"Running MeuAnalysis Constructor..."<<endl;
00075 
00076   fChain=0;
00077   fChainBD=0;
00078   fev=0;//atnu one
00079   fOutFile=0;
00080   fRec=0;//sr
00081   fRecBD=0;//sr beam data
00082   fUseAtNu=false;
00083   
00084   fEntries=0;
00085   fEntriesBD=0;
00086   fRun=100;//default
00087   fSubrun=0;
00088 
00089   if (fLoadTrees) {
00090     MSG("MeuAnalysis",Msg::kInfo)<<"Loading trees"<<endl;
00091     this->MakeChain();
00092     this->SetRootFileObjects();
00093     //this->InitialiseHitInfoVariables();
00094   
00095     //get basic run number stuff
00096     if (fRec && fChain){//sr
00097       //get the first event
00098       fChain->GetEvent(0); 
00099       NtpStRecord& ntp=(*fRec);
00100       const RecCandHeader& rec=ntp.GetHeader();
00101       MSG("MeuAnalysis",Msg::kInfo)
00102         <<"Found: run="<<rec.GetRun()
00103         <<", subrun="<<rec.GetSubRun()<<endl;
00104       fRun=rec.GetRun();
00105       fSubrun=rec.GetSubRun();
00106     }
00107     else if (fev && fChain){//atnu
00108       //get the first event
00109       fChain->GetEvent(0); 
00110       MSG("MeuAnalysis",Msg::kInfo)
00111         <<"Found: run="<<fev->Run
00112         <<", subrun="<<fev->SubRun<<endl;
00113       fRun=fev->Run;
00114       fSubrun=fev->SubRun;
00115     }
00116     else cout<<"Ahhh, no root files"<<endl;
00117   
00118   }
00119   else  {
00120     MSG("MeuAnalysis",Msg::kInfo)
00121       <<"Not loading trees in constructor"<<endl;
00122   }
00123   
00124   //setup some root stuff
00125   //include the under and overflow counts
00126   gStyle->SetOptStat(1111111);
00127   gStyle->SetOptFit(1111);
00128   gStyle->SetPalette(1);
00129   
00130   MSG("MeuAnalysis",Msg::kInfo)
00131     <<"Finished MeuAnalysis Constructor"<<endl;
00132 }

MeuAnalysis::~MeuAnalysis (  ) 

Definition at line 136 of file MeuAnalysis.cxx.

References fOutFile, Msg::kDebug, Msg::kInfo, and MSG.

00137 {
00138   MSG("MeuAnalysis",Msg::kInfo)
00139     <<"Running MeuAnalysis Destructor..."<<endl;
00140   
00141   if (fOutFile){
00142     //this makes histos disappear from the canvases 
00143     //so do it at the very end
00144     fOutFile->Close();
00145   }
00146   
00147   MSG("MeuAnalysis",Msg::kDebug)
00148     <<"Finished MeuAnalysis Destructor"<<endl;
00149 }


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, fEntries, NtpSRTrack::fidend, FilterLI(), Form(), fOutFile, fRec, InitialiseLoopVariables(), Msg::kDebug, Msg::kInfo, MAXMSG, MSG, NtpSREventSummary::ndigit, NtpSRTrack::nstrip, NtpSREventSummary::ntrack, OpenFile(), NtpSRStrip::ph1, MeuHitInfo::Plane, NtpSRStrip::plane, NtpSRVertex::plane, NtpSRStrip::planeview, NtpSRTrack::range, MeuSummary::RFid, SetLoopVariables(), NtpSRPulseHeight::sigcor, MeuHitInfo::SigCor2, NtpStRecord::stp, NtpSRTrack::stp, NtpSRTrack::stpx, NtpSRTrack::stpy, NtpSRTrack::stpz, NtpSRStrip::strip, MeuHitInfo::Strip, MeuHitInfo::TPos, NtpSRStrip::tpos, NtpStRecord::trk, MeuHitInfo::View, NtpSRTrack::vtx, NtpSRVertex::x, MeuHitInfo::X, NtpSRVertex::y, MeuHitInfo::Y, NtpSRStrip::z, NtpSRVertex::z, and MeuHitInfo::Z.

01608 {
01609   MSG("MeuAnalysis",Msg::kInfo) 
01610     <<" ** Running BasicPlots method... **"<<endl;
01611 
01612   //open the output file for the histograms
01613   fOutFile=this->OpenFile(100,"BasicPlots");
01614 
01615 
01616   TH1F* hNTrack=new TH1F("hNTrack","hNTrack",10,-1,9);
01617   hNTrack->SetTitle("NTrack");
01618   hNTrack->GetXaxis()->SetTitle("NTrack");
01619   hNTrack->GetXaxis()->CenterTitle();
01620   hNTrack->GetYaxis()->SetTitle("");
01621   hNTrack->GetYaxis()->CenterTitle();
01622   hNTrack->SetFillColor(0);
01623   //hNTrack->SetBit(TH1::kCanRebin);
01624 
01625   TH1F* hStPerTrk=new TH1F("hStPerTrk","hStPerTrk",201,-1,200);
01626   hStPerTrk->SetTitle("StPerTrk");
01627   hStPerTrk->GetXaxis()->SetTitle("StPerTrk");
01628   hStPerTrk->GetXaxis()->CenterTitle();
01629   hStPerTrk->GetYaxis()->SetTitle("");
01630   hStPerTrk->GetYaxis()->CenterTitle();
01631   hStPerTrk->SetFillColor(0);
01632   //hStPerTrk->SetBit(TH1::kCanRebin);
01633 
01634   TH1F* hTPos=new TH1F("hTPos","hTPos",200,-4,4);
01635   hTPos->SetTitle("TPos");
01636   hTPos->GetXaxis()->SetTitle("TPos");
01637   hTPos->GetXaxis()->CenterTitle();
01638   hTPos->GetYaxis()->SetTitle("");
01639   hTPos->GetYaxis()->CenterTitle();
01640   hTPos->SetFillColor(0);
01641   //hTPos->SetBit(TH1::kCanRebin);
01642 
01643   TH1F* hSigCor=new TH1F("hSigCor","hSigCor",3000,-1,15000);
01644   hSigCor->SetTitle("SigCor");
01645   hSigCor->GetXaxis()->SetTitle("SigCor");
01646   hSigCor->GetXaxis()->CenterTitle();
01647   hSigCor->GetYaxis()->SetTitle("");
01648   hSigCor->GetYaxis()->CenterTitle();
01649   hSigCor->SetFillColor(0);
01650   //hSigCor->SetBit(TH1::kCanRebin);
01651 
01652   TH1F* hView=new TH1F("hView","hView",6,-1,5);
01653   hView->SetTitle("View");
01654   hView->GetXaxis()->SetTitle("View");
01655   hView->GetXaxis()->CenterTitle();
01656   hView->GetYaxis()->SetTitle("");
01657   hView->GetYaxis()->CenterTitle();
01658   hView->SetFillColor(0);
01659   //hView->SetBit(TH1::kCanRebin);
01660 
01661   TH1F* hPlane=new TH1F("hPlane","hPlane",301,-1,300);
01662   hPlane->SetTitle("Plane");
01663   hPlane->GetXaxis()->SetTitle("Plane");
01664   hPlane->GetXaxis()->CenterTitle();
01665   hPlane->GetYaxis()->SetTitle("");
01666   hPlane->GetYaxis()->CenterTitle();
01667   hPlane->SetFillColor(0);
01668   //hPlane->SetBit(TH1::kCanRebin);
01669 
01670   TH1F* hStrip=new TH1F("hStrip","hStrip",101,-1,100);
01671   hStrip->SetTitle("Strip");
01672   hStrip->GetXaxis()->SetTitle("Strip");
01673   hStrip->GetXaxis()->CenterTitle();
01674   hStrip->GetYaxis()->SetTitle("");
01675   hStrip->GetYaxis()->CenterTitle();
01676   hStrip->SetFillColor(0);
01677   //hStrip->SetBit(TH1::kCanRebin);
01678 
01679   TH1F* hZ=new TH1F("hZ","hZ",100,-1,30);
01680   hZ->SetTitle("Z");
01681   hZ->GetXaxis()->SetTitle("Z");
01682   hZ->GetXaxis()->CenterTitle();
01683   hZ->GetYaxis()->SetTitle("");
01684   hZ->GetYaxis()->CenterTitle();
01685   hZ->SetFillColor(0);
01686   //hZ->SetBit(TH1::kCanRebin);
01687 
01688   TH1F* hR=new TH1F("hR","hR",100,-1,10);
01689   hR->SetTitle("R");
01690   hR->GetXaxis()->SetTitle("R");
01691   hR->GetXaxis()->CenterTitle();
01692   hR->GetYaxis()->SetTitle("");
01693   hR->GetYaxis()->CenterTitle();
01694   hR->SetFillColor(0);
01695   //hR->SetBit(TH1::kCanRebin);
01696 
01697   TH1F* hMaxR=new TH1F("hMaxR","hMaxR",100,-1,10);
01698   hMaxR->SetTitle("MaxR");
01699   hMaxR->GetXaxis()->SetTitle("MaxR");
01700   hMaxR->GetXaxis()->CenterTitle();
01701   hMaxR->GetYaxis()->SetTitle("");
01702   hMaxR->GetYaxis()->CenterTitle();
01703   hMaxR->SetFillColor(0);
01704   //hMaxR->SetBit(TH1::kCanRebin);
01705 
01706   TH1F* hMinR=new TH1F("hMinR","hMinR",100,-1,10);
01707   hMinR->SetTitle("MinR");
01708   hMinR->GetXaxis()->SetTitle("MinR");
01709   hMinR->GetXaxis()->CenterTitle();
01710   hMinR->GetYaxis()->SetTitle("");
01711   hMinR->GetYaxis()->CenterTitle();
01712   hMinR->SetFillColor(0);
01713   //hMinR->SetBit(TH1::kCanRebin);
01714 
01715   TH2F* hTPosVsPlane=new TH2F("hTPosVsPlane","hTPosVsPlane",
01716                               301,-1,300,100,-3,3);
01717   hTPosVsPlane->SetTitle("TPos vs Plane");
01718   hTPosVsPlane->GetXaxis()->SetTitle("Plane");
01719   hTPosVsPlane->GetXaxis()->CenterTitle();
01720   hTPosVsPlane->GetYaxis()->SetTitle("TPos");
01721   hTPosVsPlane->GetYaxis()->CenterTitle();
01722   hTPosVsPlane->SetFillColor(0);
01723   //hTPosVsPlane->SetBit(TH1::kCanRebin);
01724 
01725   TH2F* hYvsX=new TH2F("hYvsX","hYvsX",
01726                        100,-5,5,100,-5,5);
01727   hYvsX->SetTitle("Y vs X");
01728   hYvsX->GetXaxis()->SetTitle("X");
01729   hYvsX->GetXaxis()->CenterTitle();
01730   hYvsX->GetYaxis()->SetTitle("Y");
01731   hYvsX->GetYaxis()->CenterTitle();
01732   hYvsX->SetFillColor(0);
01733   //hYvsX->SetBit(TH1::kCanRebin);
01734 
01735   TH2F* hYvsXVtx=new TH2F("hYvsXVtx","hYvsXVtx",
01736                           100,-5,5,100,-5,5);
01737   hYvsXVtx->SetTitle("Y vs X");
01738   hYvsXVtx->GetXaxis()->SetTitle("X");
01739   hYvsXVtx->GetXaxis()->CenterTitle();
01740   hYvsXVtx->GetYaxis()->SetTitle("Y");
01741   hYvsXVtx->GetYaxis()->CenterTitle();
01742   hYvsXVtx->SetFillColor(0);
01743   //hYvsXVtx->SetBit(TH1::kCanRebin);
01744 
01745   TH2F* hYvsXEnd=new TH2F("hYvsXEnd","hYvsXEnd",
01746                           100,-5,5,100,-5,5);
01747   hYvsXEnd->SetTitle("Y vs X");
01748   hYvsXEnd->GetXaxis()->SetTitle("X");
01749   hYvsXEnd->GetXaxis()->CenterTitle();
01750   hYvsXEnd->GetYaxis()->SetTitle("Y");
01751   hYvsXEnd->GetYaxis()->CenterTitle();
01752   hYvsXEnd->SetFillColor(0);
01753   //hYvsXEnd->SetBit(TH1::kCanRebin);
01754 
01755   TH2F* hYvsXVtxFid=new TH2F("hYvsXVtxFid","hYvsXVtxFid",
01756                              100,-5,5,100,-5,5);
01757   hYvsXVtxFid->SetTitle("Y vs X");
01758   hYvsXVtxFid->GetXaxis()->SetTitle("X");
01759   hYvsXVtxFid->GetXaxis()->CenterTitle();
01760   hYvsXVtxFid->GetYaxis()->SetTitle("Y");
01761   hYvsXVtxFid->GetYaxis()->CenterTitle();
01762   hYvsXVtxFid->SetFillColor(0);
01763   //hYvsXVtxFid->SetBit(TH1::kCanRebin);
01764 
01765   TH2F* hYvsXEndFid=new TH2F("hYvsXEndFid","hYvsXEndFid",
01766                              100,-5,5,100,-5,5);
01767   hYvsXEndFid->SetTitle("Y vs X");
01768   hYvsXEndFid->GetXaxis()->SetTitle("X");
01769   hYvsXEndFid->GetXaxis()->CenterTitle();
01770   hYvsXEndFid->GetYaxis()->SetTitle("Y");
01771   hYvsXEndFid->GetYaxis()->CenterTitle();
01772   hYvsXEndFid->SetFillColor(0);
01773   //hYvsXEndFid->SetBit(TH1::kCanRebin);
01774 
01775   TH2F* hYvsZ=new TH2F("hYvsZ","hYvsZ",400,-1,30,100,-5,5);
01776   hYvsZ->SetTitle("Y vs Z");
01777   hYvsZ->GetXaxis()->SetTitle("Z");
01778   hYvsZ->GetXaxis()->CenterTitle();
01779   hYvsZ->GetYaxis()->SetTitle("Y");
01780   hYvsZ->GetYaxis()->CenterTitle();
01781   hYvsZ->SetFillColor(0);
01782   //hYvsZ->SetBit(TH1::kCanRebin);
01783 
01784   TH2F* hXvsZ=new TH2F("hXvsZ","hXvsZ",400,-1,30,100,-5,5);
01785   hXvsZ->SetTitle("X vs Z");
01786   hXvsZ->GetXaxis()->SetTitle("Z");
01787   hXvsZ->GetXaxis()->CenterTitle();
01788   hXvsZ->GetYaxis()->SetTitle("Y");
01789   hXvsZ->GetYaxis()->CenterTitle();
01790   hXvsZ->SetFillColor(0);
01791   //hXvsZ->SetBit(TH1::kCanRebin);
01792 
01793   TH2F* hYvsPlane=new TH2F("hYvsPlane","hYvsPlane",301,-1,300,100,-5,5);
01794   hYvsPlane->SetTitle("Y vs Z");
01795   hYvsPlane->GetXaxis()->SetTitle("Z");
01796   hYvsPlane->GetXaxis()->CenterTitle();
01797   hYvsPlane->GetYaxis()->SetTitle("Y");
01798   hYvsPlane->GetYaxis()->CenterTitle();
01799   hYvsPlane->SetFillColor(0);
01800   //hYvsPlane->SetBit(TH1::kCanRebin);
01801 
01802   TH2F* hXvsPlane=new TH2F("hXvsPlane","hXvsPlane",301,-1,300,100,-5,5);
01803   hXvsPlane->SetTitle("X vs Z");
01804   hXvsPlane->GetXaxis()->SetTitle("Z");
01805   hXvsPlane->GetXaxis()->CenterTitle();
01806   hXvsPlane->GetYaxis()->SetTitle("Y");
01807   hXvsPlane->GetYaxis()->CenterTitle();
01808   hXvsPlane->SetFillColor(0);
01809   //hXvsPlane->SetBit(TH1::kCanRebin);
01810 
01811   vector<TH1F*> vhStrips(282);
01812   vector<TH1F*> vhTPos(282);
01813   for (Int_t i=0;i<282;i++) {
01814     string sName="hStripsPlane";
01815     string sPl=Form("%d",i);
01816     sName+=sPl;
01817     //cout<<"name="<<sName<<endl;
01818     TH1F* hStrip=new TH1F(sName.c_str(),sName.c_str(),
01819                           101,-1,100);
01820     hStrip->SetTitle(sName.c_str());
01821     hStrip->GetXaxis()->SetTitle("Strip");
01822     hStrip->GetXaxis()->CenterTitle();
01823     hStrip->GetYaxis()->SetTitle("");
01824     hStrip->GetYaxis()->CenterTitle();
01825     hStrip->SetFillColor(0);
01826     //hStrip->SetBit(TH1::kCanRebin);
01827 
01828     vhStrips[i]=hStrip;
01829     //strips are 0-95(96), 4-67(64), 0-63(64)
01830       
01831     sName="hTPosPlane";
01832     sName+=sPl;
01833     TH1F* hTPos=new TH1F(sName.c_str(),sName.c_str(),200,-4,4);
01834     hTPos->SetTitle(sName.c_str());
01835     hTPos->GetXaxis()->SetTitle("TPos");
01836     hTPos->GetXaxis()->CenterTitle();
01837     hTPos->GetYaxis()->SetTitle("");
01838     hTPos->GetYaxis()->CenterTitle();
01839     hTPos->SetFillColor(0);
01840     //hTPos->SetBit(TH1::kCanRebin);
01841 
01842     vhTPos[i]=hTPos;
01843   }
01844 
01845   vector<TH1F*> vhStripsTrk(282);
01846   vector<TH1F*> vhTPosTrk(282);
01847   for (Int_t i=0;i<282;i++) {
01848     string sName="hStripsTrkPlane";
01849     string sPl=Form("%d",i);
01850     sName+=sPl;
01851     //cout<<"name="<<sName<<endl;
01852     TH1F* hStrip=new TH1F(sName.c_str(),sName.c_str(),
01853                           101,-1,100);
01854     hStrip->SetTitle(sName.c_str());
01855     hStrip->GetXaxis()->SetTitle("Strip");
01856     hStrip->GetXaxis()->CenterTitle();
01857     hStrip->GetYaxis()->SetTitle("");
01858     hStrip->GetYaxis()->CenterTitle();
01859     hStrip->SetFillColor(0);
01860     //hStrip->SetBit(TH1::kCanRebin);
01861 
01862     vhStripsTrk[i]=hStrip;
01863     //strips are 0-95(96), pl 8,10: 4-67(64), pl 7,9: 0-63(64)
01864     //plane 6 has low strips in PI part (I think 0-67)
01865     //plane 11 has high strips in PI part (I think -95)
01866 
01867     //strips TPos:
01868     //plane 6 goes -2.64 -> 1.32 m
01869     //plane 11 goes -1.32 -> 2.64 m
01870 
01871     //plane 4,8,10 goes -2.40 -> 0.24 m
01872     //plane 5,7,9 goes -0.24 -> 2.40 m
01873 
01874     //2.64 - 2.40 = 24 cm = 5.85 strips
01875     //the FI planes "stick out" by ~6 strips
01876 
01877     sName="hTPosTrkPlane";
01878     sName+=sPl;
01879     TH1F* hTPos=new TH1F(sName.c_str(),sName.c_str(),200,-4,4);
01880     hTPos->SetTitle(sName.c_str());
01881     hTPos->GetXaxis()->SetTitle("TPos");
01882     hTPos->GetXaxis()->CenterTitle();
01883     hTPos->GetYaxis()->SetTitle("");
01884     hTPos->GetYaxis()->CenterTitle();
01885     hTPos->SetFillColor(0);
01886     //hTPos->SetBit(TH1::kCanRebin);
01887 
01888     vhTPosTrk[i]=hTPos;
01889   }
01890 
01891   //reconstruction object
01892   MeuReco reco;
01893   MeuSummary s;
01894 
01898   
01899   this->InitialiseLoopVariables();  
01900   
01901   //for(Int_t entry=0;entry<10000;entry++){
01902   for(Int_t entry=0;entry<fEntries;entry++){
01903       
01904     this->SetLoopVariables(entry);
01905       
01906     NtpStRecord& ntpst=(*fRec);//Make a reference instead of a pointer
01907     
01908     //
01909     this->FilterLI(ntpst,s);
01910     continue;
01911 
01912     TClonesArray& tcaTk=(*ntpst.trk);
01913     Int_t numTrks=tcaTk.GetEntries();
01914     
01915     MAXMSG("MeuAnalysis",Msg::kInfo,10)
01916       <<"numTrks="<<numTrks
01917       <<", ndigit="<<fRec->evthdr.ndigit
01918       <<", nTrack="<<fRec->evthdr.ntrack<<endl;
01919 
01920     MeuSummary s;
01921     s.RFid=1.0;
01922     map<Int_t,MeuHitInfo> cp;
01923     
01924     hNTrack->Fill(numTrks);
01925     if (numTrks>1) continue;
01926 
01927     //Loop over tracks
01928     for (Int_t itrk=0;itrk<numTrks;itrk++){
01929       const NtpSRTrack* ptrk=
01930         dynamic_cast<NtpSRTrack*>(tcaTk[itrk]);
01931       const NtpSRTrack& trk=(*ptrk);
01932 
01933       MAXMSG("MeuAnalysis",Msg::kInfo,10)
01934         <<"range="<<trk.range<<" nstrip="<<trk.nstrip<<" ds="<<trk.ds
01935         <<" fidend.dr,dz="<<trk.fidend.dr<<" "<<trk.fidend.dz<<" end-vtx.z:"
01936         <<trk.end.z-trk.vtx.z<<" "<<trk.stpz[1]-trk.stpz[0]<<endl;
01937       
01938       hStPerTrk->Fill(trk.nstrip);
01939       hYvsXVtx->Fill(trk.vtx.x,trk.vtx.y);
01940       hYvsXEnd->Fill(trk.end.x,trk.end.y);
01941       if (trk.end.plane>5 && trk.end.plane<115){
01942         hYvsXVtxFid->Fill(trk.vtx.x,trk.vtx.y);
01943         hYvsXEndFid->Fill(trk.end.x,trk.end.y);
01944       }
01945 
01946       TClonesArray& tcaStp=(*ntpst.stp);
01947       Int_t numStps=tcaStp.GetEntries();
01948       MAXMSG("MeuAnalysis",Msg::kDebug,10)
01949         <<"numStps="<<numStps<<endl;
01950       
01951       Float_t minR=999;
01952       Float_t maxR=0;
01953 
01954       for (Int_t i=0;i<trk.nstrip;i++){
01955 
01956         //check for bug where strip index is -1
01957         if (trk.stp[i]<0) {
01958           MAXMSG("MeuAnalysis",Msg::kInfo,500)
01959             <<"Skipping strip with trk.stp[i]="<<trk.stp[i]<<endl;
01960           continue;
01961         }
01962 
01963         //const NtpSRStrip* pstp=dynamic_cast<NtpSRStrip*>(tcaStp[i]);
01964         const NtpSRStrip* pstp=
01965           dynamic_cast<NtpSRStrip*>(tcaStp[trk.stp[i]]);
01966 
01967         const NtpSRStrip& stp=(*pstp);
01968         MAXMSG("MeuAnalysis",Msg::kInfo,20)
01969           <<"Strip="<<stp.strip<<", pl="<<stp.plane
01970           <<", sigCor="<<stp.ph1.sigcor<<", tpos="<<stp.tpos<<endl;
01971         
01972         vhStripsTrk[stp.plane]->Fill(stp.strip);
01973         vhTPosTrk[stp.plane]->Fill(stp.tpos);
01974  
01975         if (stp.plane>120) continue;
01976 
01977         hTPos->Fill(stp.tpos);
01978         hSigCor->Fill(stp.ph1.sigcor);
01979         hView->Fill(stp.planeview);
01980         hPlane->Fill(stp.plane);
01981         hStrip->Fill(stp.strip);
01982         hZ->Fill(stp.z);
01983 
01984         //fill the calibpos
01985         MeuHitInfo& c=cp[stp.plane];
01986         c.Plane=stp.plane;
01987         c.View=stp.planeview;
01988         c.Strip=stp.strip;
01989         c.TPos=stp.tpos;
01990         c.SigCor2=stp.ph1.sigcor;
01991         c.X=trk.stpx[i];
01992         c.Y=trk.stpy[i];
01993         c.Z=trk.stpz[i];
01994         
01995         hTPosVsPlane->Fill(stp.plane,stp.tpos);
01996         MAXMSG("MeuAnalysis",Msg::kInfo,10)
01997           <<"x="<<trk.stpx[i]
01998           <<", y="<<trk.stpy[i]
01999           <<", z="<<trk.stpz[i]<<endl;
02000         hYvsX->Fill(trk.stpx[i],trk.stpy[i]);
02001         hYvsZ->Fill(trk.stpz[i],trk.stpy[i]);
02002         hXvsZ->Fill(trk.stpz[i],trk.stpx[i]);
02003         hYvsPlane->Fill(stp.plane,trk.stpy[i]);
02004         hXvsPlane->Fill(stp.plane,trk.stpx[i]);
02005         
02006         Float_t r=sqrt(pow((trk.stpx[i]-1.5),2)+pow(trk.stpy[i],2));
02007         if (r<minR) minR=r;
02008         if (r>maxR) maxR=r;
02009         hR->Fill(r);
02010 
02011       }
02012 
02013       //fill r histos
02014       if (maxR!=0 && minR!=999) {
02015         hMaxR->Fill(maxR);
02016         hMinR->Fill(minR);
02017       }
02018     }
02019 
02020     TClonesArray& tcaStp=(*ntpst.stp);
02021     Int_t numStps=tcaStp.GetEntries();
02022     for (Int_t i=0;i<numStps;i++){      
02023       const NtpSRStrip* pstp=dynamic_cast<NtpSRStrip*>(tcaStp[i]);
02024       const NtpSRStrip& stp=(*pstp);
02025       //MAXMSG("MeuAnalysis",Msg::kInfo,200)
02026       //  <<"Strip="<<stp.strip<<", pl="<<stp.plane<<endl;
02027         
02028       vhStrips[stp.plane]->Fill(stp.strip);
02029       vhTPos[stp.plane]->Fill(stp.tpos);
02030     }
02031 
02032   }//end of for                                       
02033   
02037 
02038   MSG("MeuAnalysis",Msg::kInfo)<<"Finished main loop"<<endl;
02039 
02040   MSG("MeuAnalysis",Msg::kInfo) 
02041     <<" ** Finished BasicPlots method **"<<endl;
02042 }

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, err(), MeuSummary::Event, NtpStRecord::evt, MeuCuts::ExtractPlInfo(), MeuSummary::FC, fEntries, MeuCuts::FillSTSumDetails(), MeuCuts::FillSTSumRecoDetails(), fRec, RecRecordImp< T >::GetHeader(), MeuSummaryWriter::GetMeuSummaryToFill(), RecDataHeader::GetRun(), RecPhysicsHeader::GetSnarl(), InitialiseLoopVariables(), Msg::kDebug, Msg::kError, Msg::kInfo, MAXMSG, MSG, NtpSRTrack::nstrip, OpenTxtFile(), MeuSummary::PC, NtpSRTrack::plane, Recalibrate(), MeuReco::Reconstruct(), MeuSummary::RFid, SetLoopVariables(), MeuSummaryWriter::SummaryTreeFill(), MeuSummaryWriter::SummaryTreeFinish(), NtpStRecord::trk, MeuSummary::VtxPlane, and MeuSummary::WinSigMap.

02927 {
02928   MSG("MeuAnalysis",Msg::kInfo) 
02929     <<" ** Running BasicReco method... **"<<endl;
02930 
02931   //open the output file for the histograms
02932   //fOutFile=this->OpenFile(100,"BasicReco");
02933 
02934   MeuSummaryWriter summaryOut;
02935 
02936   TH1F* hWinSigMap=new TH1F("hWinSigMap","hWinSigMap",3000,-1,4000);
02937   hWinSigMap->SetTitle("hWinSigMap");
02938   hWinSigMap->GetXaxis()->SetTitle("hWinSigMap");
02939   hWinSigMap->GetXaxis()->CenterTitle();
02940   hWinSigMap->SetFillColor(0);
02941   hWinSigMap->SetLineWidth(3);
02942   //hWinSigMap->SetBit(TH1::kCanRebin);
02943 
02944   Int_t FCCounter=0;
02945   Int_t PCCounter=0;
02946   Int_t TGCounter=0;
02947   Int_t badWindowCounter=0;
02948   Int_t badRecoCounter=0;
02949   Int_t goodWindowCounter=0;
02950   Int_t badFidCounter=0;
02951 
02952   ofstream& eventsTxtFile=*(this->OpenTxtFile(100,"myevents"));
02953 
02954   //reconstruction object
02955   MeuReco reco;
02956   MeuCuts cuts;
02957   
02958   //variable to turn on all the useful messages if required
02959   Msg::LogLevel_t logLevel=Msg::kDebug;
02960 
02964   
02965   this->InitialiseLoopVariables();  
02966   
02967   //for(Int_t entry=0;entry<10000;entry++){
02968   for(Int_t entry=0;entry<fEntries;entry++){
02969       
02970     this->SetLoopVariables(entry);
02971     
02972     MSG("MeuAnalysis",logLevel)
02973       <<"Entry="<<entry<<endl;
02974 
02975     //make a reference instead of a pointer
02976     NtpStRecord& ntp=(*fRec);
02977     const RecCandHeader& rec=ntp.GetHeader();
02978     MAXMSG("MeuAnalysis",Msg::kInfo,1)
02979       <<"First entry: run="<<rec.GetRun()
02980       <<", snarl="<<rec.GetSnarl()<<endl;
02981 
02982     //MeuSummary s;
02983     MeuSummary& s=summaryOut.GetMeuSummaryToFill();
02984     s.RFid=1.0;
02985     s.DistToEdgeFid=0.4;
02986     s.Event=entry;
02987     cuts.FillSTSumDetails(ntp,s,-1,-1);
02988     //a map to hold the info for each plane in the entry
02989     map<Int_t,MeuHitInfo> plInfo;
02990 
02991     //ensure there is only 1 track in the entry
02992     TClonesArray& trkTca=(*ntp.trk);
02993     Int_t numTrks=trkTca.GetEntries();
02994     if (numTrks!=1) continue;
02995     const NtpSRTrack* ptrk=dynamic_cast<NtpSRTrack*>(trkTca[0]);
02996     const NtpSRTrack& trk=(*ptrk);
02997     
02998     //cut on number of planes
02999     if (sqrt(pow(1.*trk.plane.beg-trk.plane.end,2))<15) continue;
03000 
03001     MSG("MeuAnalysis",logLevel)
03002       <<"Entry="<<entry<<", run="<<rec.GetRun()
03003       <<", snarl="<<rec.GetSnarl()
03004       <<", trk.nstrips="<<trk.nstrip<<endl;
03005 
03006     TClonesArray& evtTca=(*ntp.evt);
03007     const NtpSREvent* pevt=dynamic_cast<NtpSREvent*>(evtTca[0]);
03008     cuts.ExtractPlInfo(ntp,s,plInfo,*pevt);
03009     
03010     reco.CalcVtx(s,plInfo);
03011 
03012     //check if the track extends into the calorimeter
03013     if (s.VtxPlane>120 && s.EndPlane>120) continue;
03014 
03015     //check if track stops in spectrometer
03016     if (s.EndPlane>120) continue;
03017 
03018     MSG("MeuAnalysis",logLevel)<<"CalcFCPC..."<<endl;
03019     reco.CalcFCPC(s,plInfo);
03020 
03021     //recalculate the vtx of the tracks coming in from the spectrometer
03022     MSG("MeuAnalysis",logLevel)<<"CalcVtxSpecialND..."<<endl;
03023     reco.CalcVtxSpecialND(s,plInfo);
03024 
03025     //cut on number of planes with potential new vertex
03026     if (sqrt(pow(1.*s.VtxPlane-s.EndPlane,2))<15) continue;
03027 
03028     Bool_t TG=!s.FC && !s.PC;
03029     
03030     if (s.PC) PCCounter++;
03031     else if (s.FC) FCCounter++;
03032     else if (TG) TGCounter++;
03033     if (s.FC && s.PC) MSG("MeuAnalysis",Msg::kError)
03034       <<"Both FC and PC!"<<endl;
03035     
03036     if (s.PC) {
03037       MSG("MeuAnalysis",logLevel)<<"Found PC event..."<<endl;
03038       Bool_t goodReco=reco.Reconstruct(s,plInfo);
03039       if (goodReco){
03040         MSG("MeuAnalysis",logLevel)<<"Recalibrating..."<<endl;
03041         TClonesArray& evtTca=(*ntp.evt);
03042         const NtpSREvent* pevt=dynamic_cast<NtpSREvent*>(evtTca[0]);
03043         this->Recalibrate(ntp,s,plInfo,*pevt);
03044           
03045         MSG("MeuAnalysis",logLevel)<<"CalcWindow..."<<endl;
03046         reco.CalcWindow(s,plInfo);
03047           
03048         MSG("MeuAnalysis",logLevel)<<"CheckContainment..."<<endl;
03049         //fill the tree once event has been analysed
03050         if (s.WinSigMap<=0) {
03051           badWindowCounter++;
03052         }
03053         else if (!reco.CheckWinContainment(s,plInfo)){
03054           badFidCounter++;
03055           //MSG("MeuAnalysis",Msg::kInfo)<<"Outside fid:"<<endl;
03056         }
03057         else {
03058           goodWindowCounter++;
03059               
03060           //Bool_t gapsOk=this->CheckTrackGaps(plInfo); 
03061           //if (!gapsOk) {
03062           //  MSG("MeuAnalysis",Msg::kInfo)<<"Large gap:"<<endl;
03063           //  reco.PrintMeuHitInfo(plInfo);
03064           //}
03065               
03066           hWinSigMap->Fill(s.WinSigMap);
03067               
03068           eventsTxtFile<<entry<<endl;
03069           MAXMSG("MeuAnalysis",Msg::kInfo,10)
03070             <<"Good window: run="
03071             <<rec.GetRun()<<", snarl="<<rec.GetSnarl()
03072             <<endl;
03073               
03074           cuts.FillSTSumRecoDetails(s,plInfo);
03075           summaryOut.SummaryTreeFill(plInfo);
03076         }
03077       }
03078       else {
03079         badRecoCounter++;
03080         MSG("MeuAnalysis",logLevel)<<"Bad reco..."<<endl;
03081       }
03082     }
03083     
03084   }//end of for                                       
03085   
03089 
03090   MSG("MeuAnalysis",Msg::kInfo)<<"Finished main loop"<<endl;
03091 
03092   Double_t quantile=0.5;//quantile to compute
03093   Double_t meu=-1;
03094   hWinSigMap->GetQuantiles(1,&meu,&quantile);
03095   Double_t err=hWinSigMap->GetRMS()/sqrt(hWinSigMap->GetEntries());
03096   Double_t relErr=err/meu;
03097 
03098   MSG("MeuAnalysis",Msg::kInfo) 
03099     <<endl<<"*** Run Summary ***"<<endl
03100     <<"Raw MEU="<<meu<<" +/- "<<err<<" ("<<100.*relErr<<"%)"<<endl
03101     //<<endl
03102 
03103     <<"Total Entries="<<fEntries<<endl
03104     /*<<"Not LI="<<notLICounter<<endl
03105       <<"Not MM="<<notMMCounter<<endl
03106       <<"Good track="<<goodTrackCounter<<endl
03107       <<"Good length="<<goodLengthCounter<<endl
03108       <<"Good end & vtx="<<goodEndAndVtxCounter<<endl
03109       <<"Good XY reco="<<goodXYCounter<<endl
03110       <<endl
03111       <<"Good containmentPC="<<containmentCounter<<endl
03112       <<"  both SM hit (when PC)="<<bothSMCounter
03113       <<"   ("<<bothSM<<"% of PC)"<<endl
03114       <<"  coil hit (when PC)   ="<<coilHitCounter
03115       <<"   ("<<coil_PC<<"% of PC)"<<endl
03116       <<"  close LI (when PC)   ="<<closeLICounter
03117       <<"   ("<<closeLI<<"% of PC)"<<endl
03118       <<"  dead Chips (when PC) ="<<deadChipsCounter
03119       <<"   ("<<deadChips<<"% of PC)"<<endl
03120       <<"  boundary TF (when PC)="<<boundaryTFCounter
03121       <<"   ("<<boundaryTF<<"% of PC)"<<endl
03122       <<"  bad views (when PC)  ="<<badViewsCounter
03123       <<"   ("<<badViews<<"% of PC)"<<endl
03124       <<"  bad trk end (when PC)="<<badTrkEndCounter
03125       <<"   ("<<badTrkEnd<<"% of PC)"<<endl
03126       <<"  bad trk gap (when PC)="<<badTrkGapsCounter
03127       <<"   ("<<badTrkGaps<<"% of PC)"<<endl
03128       <<"  bad traceZ (when PC) ="<<badTraceZCounter
03129       <<"   ("<<badTraceZ<<"% of PC)"<<endl
03130       <<"  bad window (when PC) ="<<badWindowCounter
03131       <<"   ("<<badWin<<"% of PC)"<<endl
03132       <<"Good containmentFC="<<containmentCounterFC<<endl
03133       <<"  both SM hit (when FC)="<<bothSMCounterFC
03134       <<"   ("<<bothSMFC<<"% of FC)"<<endl
03135       <<"  coil hit (when FC)   ="<<coilHitCounterFC
03136       <<"   ("<<coil_FC<<"% of FC)"<<endl
03137       <<"Good containmentTG="<<containmentCounterTG<<endl
03138       <<"  both SM hit (when TG)="<<bothSMCounterTG
03139       <<"   ("<<bothSMTG<<"% of TG)"<<endl
03140       <<"  coil hit (when TG)   ="<<coilHitCounterTG
03141       <<"   ("<<coil_TG<<"% of TG)"<<endl
03142       <<endl
03143       <<"Bad charge pos moment="<<momentCounter<<endl
03144       <<"Too high StripsPerPlane="<<stPerPlCounter<<endl
03145       <<"  Either of above two="<<stPlOrMomentCounter<<endl
03146       <<endl
03147       <<"MC Info:"<<endl
03148       <<"  Muon lowest energy >0.6 GeV="<<lowEnCounter<<endl
03149       <<"  Muon lowest energy >5.0 GeV="<<lowEn5Counter<<endl
03150       <<endl
03151     */
03152     <<"Containment summary:"<<endl
03153     <<"  PC="<<PCCounter<<endl
03154 //    <<"     (PC - all cuts="
03156     //closeLICounter-deadChipsCounter-boundaryTFCounter-
03157     //badViewsCounter-badTrkEndCounter-badTrkGapsCounter-badTraceZCounter-
03158     //badWindowCounter<<")"<<endl
03159     <<"  FC="<<FCCounter<<endl//<<"     (FC-coil-bothSM="
03160     //<<FCCounter-coilHitCounterFC-bothSMCounterFC<<")"<<endl
03161     <<"  TG="<<TGCounter<<endl//<<"     (TG-coil-bothSM="
03162 //    <<TGCounter-coilHitCounterTG-bothSMCounterTG<<")"<<endl
03163     <<"  FC+PC+TG="<<FCCounter+PCCounter+TGCounter<<endl
03164     <<endl
03165     <<"Total PC with good track window = "<<goodWindowCounter<<endl
03166     <<"              bad track window  = "<<badWindowCounter<<endl
03167     <<"              bad fiducial vol  = "<<badFidCounter<<endl
03168     <<"              bad reco          = "<<badRecoCounter<<endl
03169     <<endl;
03170 
03171   //finish the summary tree and output it
03172   summaryOut.SummaryTreeFinish();
03173 
03174   if (fRec) delete fRec;
03175 
03176   MSG("MeuAnalysis",Msg::kInfo) 
03177     <<" ** Finished BasicReco method **"<<endl;
03178 }

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

Definition at line 2087 of file MeuAnalysis.cxx.

References Msg::kInfo, MAXMSG, NtpSRTrack::nstrip, NtpSRStrip::plane, NtpStRecord::stp, NtpSRTrack::stp, and NtpStRecord::trk.

02088 {
02089   TClonesArray& tcaTk=(*ntp.trk);
02090   Int_t numTrks=tcaTk.GetEntries();
02091   
02092   //Loop over tracks
02093   for (Int_t itrk=0;itrk<numTrks;itrk++){
02094     const NtpSRTrack* ptrk=
02095       dynamic_cast<NtpSRTrack*>(tcaTk[itrk]);
02096     const NtpSRTrack& trk=(*ptrk);
02097       
02098     TClonesArray& tcaStp=(*ntp.stp);
02099       
02100     //this loop is just over the tracked strips
02101     for (Int_t i=0;i<trk.nstrip;i++){
02102         
02103       //check for bug where strip index is -1
02104       if (trk.stp[i]<0) {
02105         MAXMSG("MeuAnalysis",Msg::kInfo,500)
02106           <<"Skipping strip with trk.stp[i]="<<trk.stp[i]<<endl;
02107         continue;
02108       }
02109 
02110       const NtpSRStrip* pstp=
02111         dynamic_cast<NtpSRStrip*>(tcaStp[trk.stp[i]]);
02112       const NtpSRStrip& stp=(*pstp);
02113       //MAXMSG("MeuAnalysis",Msg::kInfo,200)
02114       //  <<"Strip="<<stp.strip<<", pl="<<stp.plane
02115       //  <<", sigCor="<<stp.ph1.sigcor<<", tpos="<<stp.tpos<<endl;
02116 
02117       Int_t pl=stp.plane;
02118 
02119       if (pl>120) continue;
02120     }
02121   }
02122 }

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

Definition at line 4170 of file MeuAnalysis.cxx.

References it, Msg::kDebug, Msg::kVerbose, MAXMSG, AtmosStrip::Plane, AtmosStrip::QPEcorr, and AtmosEvent::StripList.

04171 {
04172   Float_t sumSigCorEast=0;
04173   Float_t sumSigCorWest=0;
04174   
04175   map<Int_t,Int_t> hpp;
04176   
04177   TClonesArray& strips=(*ev.StripList);
04178   Int_t numStrips=strips.GetEntries();
04179   if (numStrips<=0) return false;//pass the event, not LI
04181   //loop over the strips in the snarl
04183   for (Int_t hit=0;hit<numStrips;hit++){
04184     const AtmosStrip* strip=dynamic_cast<AtmosStrip*>(strips[hit]);
04185     const AtmosStrip& st=(*strip);
04186     
04187     sumSigCorEast+=st.QPEcorr[0];
04188     sumSigCorWest+=st.QPEcorr[1];
04189     
04190     hpp[st.Plane]++;
04191   }//end of loop over strips
04192   
04193   Float_t avHpp=0;
04194   vector<Float_t> pb(8);//to store fractions hit in each pb region
04195   
04196   //loop over all the planes flashed
04197   for (map<Int_t,Int_t>::iterator hppIt=hpp.begin();
04198        hppIt!=hpp.end();++hppIt){
04199     avHpp+=hppIt->second;
04200     
04201     Int_t pl=hppIt->first;
04202     
04203     //work out the fraction of the planes flashed by each pb
04204     //are actually hit
04205     if (pl>=1 && pl<=64) pb[0]+=1./64;//0=bookend
04206     else if (pl>=65 && pl<=128) pb[1]+=1./64;
04207     else if (pl>=129 && pl<=192) pb[2]+=1./64;
04208     else if (pl>=193 && pl<=248) pb[3]+=1./56;//249=bookend
04209     else if (pl>=250 && pl<=313) pb[4]+=1./64;
04210     else if (pl>=314 && pl<=377) pb[5]+=1./64;
04211     else if (pl>=378 && pl<=441) pb[6]+=1./64;
04212     else if (pl>=442 && pl<=485) pb[7]+=1./44;
04213   }
04214   if (hpp.size()>0) avHpp/=hpp.size();
04215   Float_t sumSig=sumSigCorEast+sumSigCorWest;
04216   Float_t minSigCorEW = TMath::Min(sumSigCorEast, sumSigCorWest);
04217   Float_t asym=0;
04218   if (sumSig) asym = (sumSig - 2 * minSigCorEW) / sumSig;
04219 
04220   Float_t fractFlashed1=0;
04221   Float_t fractFlashed2=0;
04222   
04223   for (vector<Float_t>::iterator it=pb.begin();
04224        it!=pb.end();++it){
04225     Float_t fract=*it;
04226     if (fract>fractFlashed1){
04227       //copy the 1 into 2
04228       fractFlashed2=fractFlashed1;
04229       //then set new highest fractFlashed
04230       fractFlashed1=fract;
04231     }
04232     MAXMSG("MeuAnalysis",Msg::kVerbose,200) 
04233       <<"fract="<<fract<<", f1="<<fractFlashed1
04234       <<", f2="<<fractFlashed2<<endl;
04235   }
04236   
04237   Float_t ratioFlashed=0;
04238   if (fractFlashed1>0) ratioFlashed=fractFlashed2/fractFlashed1;
04239   
04240   Bool_t li=false;
04241   
04242   //check if LI
04243   if ( (asym>0.5 || minSigCorEW>1.7e6) && 
04244        avHpp>3 && ratioFlashed<0.05 &&
04245        fractFlashed1>0.85 && fractFlashed2<0.1 ){
04246     li=true;
04247     
04248     MAXMSG("MeuAnalysis",Msg::kDebug,200) 
04249       <<"LI Event:"<<endl
04250       <<"  Asymmetry="<<asym<<endl
04251       <<"  Av. HPP  ="<<avHpp<<endl
04252       <<"  Fract Fl1="<<fractFlashed1<<", Fract Fl2="<<fractFlashed2
04253       <<"  Ratio Fla="<<ratioFlashed<<endl;
04254   }
04255   
04256   return li;
04257 }

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

Definition at line 4025 of file MeuAnalysis.cxx.

References MeuSummary::Detector, MuELoss::e, it, Msg::kDebug, Detector::kFar, Msg::kInfo, MAXMSG, NtpSRStrip::ph0, NtpSRStrip::ph1, NtpSRStrip::plane, NtpSRPulseHeight::sigcor, NtpStRecord::stp, NtpSRStrip::time0, and NtpSRStrip::time1.

Referenced by BasicPlots(), and MakeSummaryTreeWithAtNu().

04027 {
04028   //only analyse if FD
04029   if (s.Detector!=Detector::kFar) return false;
04030 
04031   const Bool_t cleanTheEvent=false;
04032 
04033   Float_t sumSigCorEast=0;
04034   Float_t sumSigCorWest=0;
04035 
04036   map<Int_t,Int_t> hpp;  // <plane, # of hits>
04037   multiset<Double_t> times;
04038   
04039   Int_t numStrips=ntp.stp->GetEntries();
04040   if (numStrips<=0) return false;//pass the event, not LI
04042   //loop over the strips in the snarl
04044   for (Int_t hit=0;hit<numStrips;hit++){
04045     NtpSRStrip* strip = dynamic_cast<NtpSRStrip*>((*ntp.stp)[hit]);
04046 
04047     Double_t time=strip->time0;
04048     if (time<0 || time>1) time=strip->time1;
04049 
04050     if (time>0 && time<=1) {
04051       times.insert(time);
04052     }
04053     //else just don't put the time in the map - clearly rubbish
04054   }//end of loop over strips
04055 
04056   //get the median time from the map
04057   multiset<Double_t>::const_iterator it=times.begin();
04058   advance(it,times.size()/2);
04059   Double_t medianTime=*it;
04060   MAXMSG("MeuAnalysis",Msg::kDebug,100)
04061     <<"FilterLI: Median time="<<medianTime<<endl;
04062 
04063   //second loop
04064   for (Int_t hit=0;hit<numStrips;hit++){
04065     NtpSRStrip* strip = dynamic_cast<NtpSRStrip*>((*ntp.stp)[hit]);
04066 
04067     if (cleanTheEvent) {
04068       if (strip->ph0.sigcor > 0){
04069         if (strip->time0<medianTime-200e-9 || 
04070             strip->time0>medianTime+200e-9) continue;
04071       }
04072       if (strip->ph1.sigcor > 0){
04073         if (strip->time1<medianTime-200e-9 || 
04074             strip->time1>medianTime+200e-9) continue;
04075       }
04076    } // if (cleanTheEvent)
04077 
04078     sumSigCorEast+=strip->ph0.sigcor;
04079     sumSigCorWest+=strip->ph1.sigcor;
04080 
04081     hpp[strip->plane] += 1;
04082   }//end of loop over strips
04083 
04084   Float_t avHpp=0;
04085   vector<Float_t> pb(8);//to store fractions hit in each pb region
04086 
04087   //loop over all the planes flashed
04088   for (map<Int_t,Int_t>::iterator hppIt=hpp.begin();
04089        hppIt!=hpp.end();++hppIt){
04090     avHpp+=hppIt->second;
04091 
04092     Int_t pl=hppIt->first;
04093 
04094     if (pl>=1 && pl<=64) pb[0]+=1./64;//0=bookend
04095     else if (pl>=65 && pl<=128) pb[1]+=1./64;
04096     else if (pl>=129 && pl<=192) pb[2]+=1./64;
04097     else if (pl>=193 && pl<=248) pb[3]+=1./56;//249=bookend
04098     else if (pl>=250 && pl<=313) pb[4]+=1./64;
04099     else if (pl>=314 && pl<=377) pb[5]+=1./64;
04100     else if (pl>=378 && pl<=441) pb[6]+=1./64;
04101     else if (pl>=442 && pl<=485) pb[7]+=1./44;
04102 
04103     MAXMSG("MeuAnalysis",Msg::kDebug,2000)
04104       <<"plane="<<pl<<", hpp="<<hppIt->second<<endl;
04105   }
04106   if (hpp.size()>0) avHpp/=hpp.size();
04107   Float_t sumSig=sumSigCorEast+sumSigCorWest;
04108   Float_t minSigCorEW = TMath::Min(sumSigCorEast, sumSigCorWest);
04109   Float_t asym=0;
04110   
04111   if (sumSig) asym = (sumSig - 2 * minSigCorEW) / sumSig;
04112 
04113   Float_t fractFlashed1=0;
04114   Float_t fractFlashed2=0;
04115 
04116   for (vector<Float_t>::iterator it=pb.begin();
04117        it!=pb.end();++it){
04118     Float_t fract=*it;
04119     if (fract>fractFlashed1){
04120       //copy the 1 into 2
04121       fractFlashed2=fractFlashed1;
04122       //then set new highest fractFlashed
04123       fractFlashed1=fract;
04124     }
04125     //make sure to get fractFlashed2 if fractFlashed2<fract<frFla1
04126     else if (fract>fractFlashed2) {
04127       fractFlashed2=fract;
04128     }
04129 
04130     MAXMSG("MeuAnalysis",Msg::kDebug,200)
04131       <<"fract="<<fract<<", f1="<<fractFlashed1
04132       <<", f2="<<fractFlashed2<<endl;
04133   }
04134 
04135   Float_t ratioFlashed=0;
04136   if (fractFlashed1>0) ratioFlashed=fractFlashed2/fractFlashed1;
04137 
04138   Bool_t li=false;
04139 
04140   // check if LI
04141   // don't check the asymmetry if the VA is saturating
04142   if ( (asym>0.5 || minSigCorEW>1.7e6) && 
04143        avHpp>3 && ratioFlashed<0.05 &&
04144        fractFlashed1>0.85 && fractFlashed2<0.1 ){
04145     li=true;
04146     
04147     MAXMSG("MeuAnalysis",Msg::kInfo,10)
04148       <<"LI Event:"<<endl
04149       <<"  Asymmetry="<<asym<<endl
04150       <<"  Av. HPP  ="<<avHpp<<endl
04151       <<"  Fract Fl1="<<fractFlashed1
04152       <<", Fract Fl2="<<fractFlashed2
04153       <<"  Ratio Fla="<<ratioFlashed<<endl;
04154   }
04155   else{
04156     MAXMSG("MeuAnalysis",Msg::kDebug,100)
04157       <<"Non-LI Event:"<<endl
04158       <<"  Asymmetry="<<asym<<endl
04159       <<"  Av. HPP  ="<<avHpp<<endl
04160       <<"  Fract Fl1="<<fractFlashed1
04161       <<", Fract Fl2="<<fractFlashed2
04162       <<"  Ratio Fla="<<ratioFlashed<<endl;
04163   }
04164   
04165   return li;
04166 }

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

Definition at line 375 of file MeuAnalysis.cxx.

References Form(), and gSystem().

Referenced by MakeChain().

00376 {
00377   //check if there is actually a wildcard
00378   if (!TString(wildcardString.c_str()).MaybeWildcard()) {
00379     return wildcardString;
00380   }
00381   
00382   //this was ripped off from ROOT's TChain::AddFile
00383 
00384   vector<string> fileList;
00385 
00386   //wildcarding used in name
00387   TString basename(wildcardString.c_str());
00388 
00389   Int_t dotslashpos = basename.Index(".root/");
00390   TString behind_dot_root;
00391   if (dotslashpos>=0) {
00392     // Copy the tree name specification
00393     behind_dot_root = basename(dotslashpos+6,basename.Length()-dotslashpos+6);
00394     // and remove it from basename
00395     basename.Remove(dotslashpos+5);
00396   }
00397 
00398   Int_t slashpos = basename.Last('/');
00399   TString directory;
00400   if (slashpos>=0) {
00401     directory = basename(0,slashpos); // Copy the directory name
00402     basename.Remove(0,slashpos+1);      // and remove it from basename
00403   } else {
00404     directory = gSystem->WorkingDirectory();
00405   }
00406 
00407   const char *file;
00408   void *dir = gSystem->OpenDirectory(gSystem->ExpandPathName(directory.Data()));
00409 
00410   if (dir) {
00411     //create a TList to store the file names (not yet sorted)
00412     TList l;
00413     TRegexp re(basename,kTRUE);
00414     while ((file = gSystem->GetDirEntry(dir))) {
00415       if (!strcmp(file,".") || !strcmp(file,"..")) continue;
00416       TString s = file;
00417       if ( (basename!=file) && s.Index(re) == kNPOS) continue;
00418       l.Add(new TObjString(file));
00419     }
00420     gSystem->FreeDirectory(dir);
00421     //sort the files in alphanumeric order
00422     l.Sort();
00423     TIter next(&l);
00424     TObjString *obj;
00425     while ((obj = (TObjString*)next())) {
00426       file = obj->GetName();
00427       if (behind_dot_root.Length()!=0){
00428         string fileName=Form("%s/%s/%s",directory.Data(),
00429                              file,behind_dot_root.Data());
00430         fileList.push_back(fileName);
00431         //nf += AddFile(Form("%s/%s/%s",directory.Data(),file,behind_dot_root.Data()),kBigNumber);
00432       }
00433       else {
00434         string fileName=Form("%s/%s",directory.Data(),file);
00435         fileList.push_back(fileName);
00436         //nf += AddFile(Form("%s/%s",directory.Data(),file),kBigNumber);
00437       }
00438     }
00439     l.Delete();
00440   }
00441 
00442   //check if any files were found
00443   if (fileList.begin()!=fileList.end()){
00444     cout<<"Used wildcard expansion to find first file name="
00445         <<*fileList.begin()<<endl;
00446     return *fileList.begin();
00447   }
00448   else {
00449     return "";//null string
00450   }
00451 }

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

Definition at line 4261 of file MeuAnalysis.cxx.

References MeuSummary::Detector, it, Msg::kDebug, Detector::kFar, Msg::kInfo, Msg::kVerbose, Msg::kWarning, MAXMSG, MSG, VldTimeStamp::Print(), and MeuSummary::TimeSec.

Referenced by MakeSummaryTreeWithAtNu().

04262 {
04263   static Bool_t firstCall=true;
04264   static Bool_t foundFile=false;
04265   static vector <Float_t> vTemperatures;
04266   static Int_t firstTime=-1;
04267   static Int_t lastTime=-1;
04268 
04269   //only get the temperature for the FD at present
04270   if (s.Detector!=Detector::kFar) {
04271     MAXMSG("MeuAnalysis",Msg::kInfo,1)
04272       <<"Not retrieving the temperature for dt="<<s.Detector<<endl;
04273     return 0;
04274   }
04275   
04276   if (firstCall){
04277     char* envVariable=getenv("MEUDATA");
04278     if (envVariable==NULL){
04279       MAXMSG("MeuAnalysis",Msg::kWarning,1)
04280         <<"Can't find ENV variable \"MEUDATA\" to"
04281         <<" open temperatures file"<<endl;
04282       firstCall=false;
04283       return 0;
04284     }
04285     string sEnv=envVariable;
04286 
04287     string sTemperatureFileName=sEnv+"/temps.dat";
04288     MSG("MeuAnalysis",Msg::kInfo)
04289       <<"GetTemperature: Looking for file: "
04290       <<sTemperatureFileName<<endl;
04291     ifstream temperatureFile(sTemperatureFileName.c_str());
04292     if(!temperatureFile){
04293       MAXMSG("MeuAnalysis",Msg::kWarning,1)
04294         <<"Can't find file of temperatures"<<endl;
04295       firstCall=false;
04296       foundFile=false;
04297       return 0;
04298     }
04299 
04300     foundFile=true;
04301 
04302     map<Int_t,Float_t> temperatures;
04303 
04304     MSG("MeuAnalysis",Msg::kInfo)
04305       <<"Loading temperatures from txt file..."<<endl;
04306     if (temperatureFile){
04307       Float_t Tf=-1;
04308       Float_t Tc=-1;
04309       Int_t time=0;
04310       while(temperatureFile>>time>>Tf) {
04311         Tc=(Tf-32)/1.8;
04312         temperatures[time]=Tc;
04313       }
04314     }
04315     
04316     firstTime=temperatures.begin()->first;
04317     lastTime=(--temperatures.end())->first;
04318     MSG("MeuAnalysis",Msg::kInfo)
04319       <<"Found first time="<<firstTime<<", lastTime="<<lastTime<<endl;
04320     VldTimeStamp first(firstTime,0);
04321     first.Print();
04322     VldTimeStamp last(lastTime,0);
04323     last.Print();
04324     
04325     Int_t lastIndex=(lastTime-firstTime)/300;
04326     MSG("MeuAnalysis",Msg::kInfo)<<"Size of table="<<lastIndex<<endl;
04327     vTemperatures.reserve(lastIndex+1);
04328     for (Int_t i=0;i<lastIndex+1;i++) vTemperatures.push_back(0);
04329     
04330     //loop and fill vector
04331     for (map<Int_t,Float_t>::iterator it=temperatures.begin();
04332          it!=temperatures.end();++it){
04333       Int_t time=it->first;
04334       Int_t index=(time-firstTime)/300;//every 5 mins
04335       MAXMSG("MeuAnalysis",Msg::kDebug,50)
04336         <<"time="<<it->first
04337         <<", mins since first="<<(time-firstTime)/300
04338         <<", temperature="<<it->second<<endl;
04339       vTemperatures[index]=it->second;
04340     } 
04341   }
04342   firstCall=false;
04343   
04344   if (!foundFile) return 0;
04345 
04346   Int_t ind=(s.TimeSec-firstTime)/300;
04347   Float_t Tc=0;
04348   Int_t lastIndex=(lastTime-firstTime)/300;
04349   Int_t timeSkipped=0;
04350 
04351   if (s.TimeSec>=firstTime && s.TimeSec<=lastTime){
04352     while (ind<=lastIndex && ind>0){
04353       MAXMSG("MeuAnalysis",Msg::kVerbose,500)
04354         <<"500 mins T="<<vTemperatures[ind]<<", ind="<<ind<<endl;
04355       if (vTemperatures[ind]>0){
04356         MAXMSG("MeuAnalysis",Msg::kDebug,50)
04357           <<"Found good T="<<vTemperatures[ind]<<endl;
04358         Tc=vTemperatures[ind];
04359         break;
04360       }
04361       timeSkipped++;
04362       ind++;
04363     }
04364   }
04365   else MSG("MeuAnalysis",Msg::kWarning)<<"Time is out of range"<<endl;
04366 
04367   //check that not too big a time (1 hour) was skipped
04368   if (timeSkipped<=12) return Tc;
04369   else {
04370     MAXMSG("MeuAnalysis",Msg::kInfo,50)
04371       <<"Too big a gap (>1 hour) in temperature reading"
04372       <<" to be reliable ("<<timeSkipped<<"*5 mins)."
04373       <<" Returning 0 degC"<<endl;
04374     return 0;
04375   }
04376 }

void MeuAnalysis::InitialiseLoopVariables (  )  [private]

Definition at line 762 of file MeuAnalysis.cxx.

References Msg::kDebug, and MSG.

Referenced by BasicPlots(), BasicReco(), MakeSummaryTreeWithAtNu(), MakeSummaryTreeWithNtpSt(), N_1Plots(), SnarlList(), SpillPlots(), and TestEventLoop().

00763 {
00764   MSG("MeuAnalysis",Msg::kDebug)<<"Initialising loop variables..."<<endl;
00765 
00766   MSG("MeuAnalysis",Msg::kDebug)<<"Initialisation complete"<<endl;
00767 }

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

Definition at line 164 of file MeuAnalysis.cxx.

References fLoadTrees, Msg::kInfo, and MSG.

00165 {
00166   MSG("MeuAnalysis",Msg::kInfo)
00167     <<"Setting fLoadTrees="<<loadTrees<<endl;
00168   fLoadTrees=loadTrees;
00169 }

void MeuAnalysis::MakeChain (  )  [private]

Definition at line 455 of file MeuAnalysis.cxx.

References exit(), fChain, fChainBD, fUseAtNu, GetFirstFileName(), Msg::kFatal, Msg::kInfo, MakeFileList(), MSG, and ObjectExistsInFile().

00456 {
00457   //get the files to open
00458   vector<string> fileList=this->MakeFileList();
00459 
00460   TFile* firstFile=0;
00461   // case with one single file
00462   if (!TString(*fileList.begin()).MaybeWildcard()) {
00463     firstFile=TFile::Open((*fileList.begin()).c_str(),"READ");
00464   }
00465   else {
00466     MSG("MeuAnalysis",Msg::kInfo)
00467       <<"Found a wildcard present in list of file names: "
00468       <<*fileList.begin()
00469       <<endl<<"Looking for first filename using wildcard..."<<endl;
00470     string fileName=this->GetFirstFileName(*fileList.begin());
00471     if (fileName!="") {
00472       firstFile=TFile::Open(fileName.c_str(),"READ");
00473     }
00474     else {
00475       MSG("MeuAnalysis",Msg::kInfo)
00476         <<"No files found matching:"<<endl
00477         <<*fileList.begin()<<", will exit here..."<<endl;
00478       exit(1);
00479     }
00480   }
00481   
00482   //create a chain
00483   if (false && fUseAtNu) {
00484     if (this->ObjectExistsInFile(firstFile,"ntp")) {
00485       fChain=new TChain("ntp");
00486     }
00487   }
00488   else {
00489     //check if the AtmosEvent ntuple exists in the file
00490     if (this->ObjectExistsInFile(firstFile,"ntp")) {
00491       fChain=new TChain("ntp");
00492       fUseAtNu=true;//set this to true
00493     }
00494     //else check if the SR ntuple exists
00495     else if (this->ObjectExistsInFile(firstFile,"NtpSt")) {
00496       fChain=new TChain("NtpSt");
00497     }
00498     else {
00499       MSG("MeuAnalysis",Msg::kInfo)
00500         <<"Not creating NtpSt TChain because it does not exist in file"
00501         <<endl;
00502     }
00503 
00504     //check in BD ntuple exists
00505     if (this->ObjectExistsInFile(firstFile,"NtpBDLite")) {
00506       fChainBD=new TChain("NtpBDLite");//have to load libraries too
00507     }
00508     else {
00509       MSG("MeuAnalysis",Msg::kInfo)
00510         <<"Not creating NtpBDLite TChain because it does not exist"
00511         <<" in file"<<endl;
00512     }
00513   }
00514   
00515   firstFile->Close();
00516   if (firstFile) delete firstFile;
00517   firstFile=0;
00518 
00519   Int_t nf=0;
00520   Int_t nfBD=0;
00521   //add the files to the chain
00522   for (vector<string>::iterator file=fileList.begin();
00523        file!=fileList.end();++file){
00524     
00525     //test if file already exists
00526     ifstream openOk((*file).c_str()); 
00527 
00528     //check if a wildcard was used because ifstream can't open wildcards
00529     Int_t stars=(*file).find("*");
00530     Int_t quest=(*file).find("?");
00531 
00532     //check if file existed
00533     if (!openOk && !(quest>=0 || stars>=0)){
00534       MSG("MeuAnalysis",Msg::kInfo)
00535         <<endl<<endl
00536         <<"***********************************************************"
00537         <<endl<<"Can't find file="<<*file<<endl
00538         <<"Note: you can't use '~/'. It has to be the full path"<<endl
00539         <<"***********************************************************"
00540         <<endl<<endl
00541         <<"Exiting here!"<<endl;
00542       exit(0);
00543     }
00544     
00545     MSG("MeuAnalysis",Msg::kInfo)<<"Adding file(s)="<<*file<<endl;
00546     nf+=fChain->Add((*file).c_str());
00547     if (fChainBD) nfBD+=fChainBD->Add((*file).c_str());
00548     //fChain->Print();//very verbose
00549   }
00550 
00551   if(nf==0){
00552     MSG("MeuAnalysis",Msg::kFatal)
00553       <<endl<<endl
00554       <<"*************************************************************"
00555       <<endl<<"No *.root files found"<<endl
00556       <<"Please set MEUDATA to the directory containing the"
00557       <<" *.root files"<<endl
00558       <<"Or give the txt file containing the files to be input"<<endl
00559       <<"Note: If more than one file is found they will be"
00560       <<" concatenated in a TChain and treated as one"<<endl
00561       <<"*************************************************************"
00562       <<endl<<endl<<"Program will exit here"<<endl;
00563     exit(0);
00564   }
00565   
00566   MSG("MeuAnalysis",Msg::kInfo)
00567     <<endl<<"Added "<<nf<<" file(s) to TChain"<<endl;
00568   MSG("MeuAnalysis",Msg::kInfo)
00569     <<endl<<"Added "<<nfBD<<" beam data (BD) file(s) to TChain"
00570     <<endl;
00571   
00572   if (fChain) {
00573     MSG("MeuAnalysis",Msg::kInfo) 
00574       <<"Ntuple information:"<<endl;
00575     fChain->Show(0);
00576     
00577     if (fChainBD) {
00578       MSG("MeuAnalysis",Msg::kInfo) 
00579         <<"NtpBDLite information:"<<endl;
00580       fChainBD->Show(0);
00581     }
00582   }
00583   
00584   MSG("MeuAnalysis",Msg::kInfo)
00585     <<endl<<"Added "<<nfBD<<" beam data (BD) file(s) to TChain."
00586     <<endl;
00587   MSG("MeuAnalysis",Msg::kInfo)
00588     <<endl<<"Analysing "<<nf<<" file(s). Reading from disk..."<<endl;
00589 }

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

Check the fInputFileName first then check the env variable

Definition at line 259 of file MeuAnalysis.cxx.

References exit(), fInputFileName, Msg::kDebug, Msg::kFatal, Msg::kInfo, and MSG.

Referenced by MakeChain().

00260 {
00262 
00263   vector<string> fileList;
00264   
00265   if (fInputFileName!=""){
00266 
00267     Int_t findGives=fInputFileName.find(".root");
00268     MSG("MeuAnalysis",Msg::kInfo)
00269       <<"Checking input string for \"*.root\", position in string="
00270       <<findGives<<endl;
00271 
00272     if (findGives>0){
00273       //add the file direct to the list
00274       fileList.push_back(fInputFileName);
00275     }
00276     else{//treat the file as a list of root files 
00277       ifstream inputFile(fInputFileName.c_str());
00278       
00279       //check if file exists
00280       if (inputFile){
00281         //variables to hold input from file
00282         string file="";
00283         
00284         //read in from the text file and fill objects
00285         while(inputFile>>file) {
00286           MSG("MeuAnalysis",Msg::kDebug)
00287             <<"Found input file name="<<file<<endl;
00288           fileList.push_back(file);
00289         }
00290         MSG("MeuAnalysis",Msg::kDebug)
00291           <<"Files names found in txt file="<<fileList.size()<<endl;
00292       }
00293       else{
00294         MSG("MeuAnalysis",Msg::kFatal)
00295           <<endl<<endl
00296           <<"**********************************************************"
00297           <<endl<<"Input txt file of file names does not exist!"<<endl
00298           <<"InputFileName="<<fInputFileName<<endl
00299           <<"**********************************************************"
00300           <<endl<<endl<<"Program will exit here"<<endl;
00301         exit(0);
00302       }
00303     }
00304   }
00305   //return the fileList if files were found
00306   if (fileList.size()>0) return fileList;
00307   
00308   //Check the env variable to find files
00309   char* envVariable=getenv("MEUDATA");
00310   if (envVariable==NULL){
00311     MSG("MeuAnalysis",Msg::kFatal)
00312       <<endl<<endl
00313       <<"*************************************************************"
00314       <<endl<<"Environmental variable MEUDATA not set!"<<endl
00315       <<"Please set MEUDATA to the directory containing the"
00316       <<" ntuple root files"<<endl
00317       <<"Note: If more than one file is found they will be"
00318       <<" concatenated and treated as one"<<endl
00319       <<"*************************************************************"
00320       <<endl<<endl<<"Program will exit here"<<endl;
00321     exit(0);
00322   }
00323 
00324   string sEnv="";
00325   sEnv=envVariable;
00326   sEnv+="/*.root";
00327   MSG("MeuAnalysis",Msg::kInfo)
00328     <<"Looking for *.root files using the env variable"<<endl
00329     <<"MEUDATA="<<sEnv<<endl;    
00330   fileList.push_back(sEnv);
00331 
00332   return fileList;
00333 }

void MeuAnalysis::MakeSummaryTree (  ) 

Definition at line 4380 of file MeuAnalysis.cxx.

References fUseAtNu, MakeSummaryTreeWithAtNu(), and MakeSummaryTreeWithNtpSt().

04381 {
04382   //this method just adds a level of indirection so that you only 
04383   //need a single macro to make the summary tree
04384 
04385   if (fUseAtNu){
04386     this->MakeSummaryTreeWithAtNu();
04387   }
04388   else{
04389     this->MakeSummaryTreeWithNtpSt();
04390   }
04391 }

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, fEntries, MeuCuts::FillEnVsDist(), MeuCuts::FillEnVsDist2(), MeuHistos::FillMeuHistos(), MeuCuts::FillSTSumDetails(), MeuCuts::FillSTSumRecoDetails(), MeuCuts::FilterBadDistEndStrip(), MeuCuts::FilterBadEndGaps(), MeuCuts::FilterBadTrackEnd(), MeuCuts::FilterBadXY(), FilterLI(), fRun, fSubrun, MeuSummaryWriter::GetMeuSummaryToFill(), GetTemperature(), Munits::GeV, MeuCounter::goodContainmentCounter, MeuCounter::goodWindowCounter, MeuCounter::goodXYCounter, InitialiseLoopVariables(), it, Msg::kDebug, Msg::kError, Msg::kInfo, MeuCounter::longTrackCounter, MeuHitInfo::LPos, Munits::m, 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, 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().

04396 {
04397   MSG("MeuAnalysis",Msg::kInfo) 
04398     <<" ** Running MakeSummaryTreeWithAtNu method... **"<<endl;
04399 
04400   //open the output file for the histograms
04401   //fOutFile=this->OpenFile(100,"MakeSummaryTreeFar");
04402 
04403   MeuSummaryWriter summaryOut;
04404   summaryOut.SummaryTreeSetup(fRun,fSubrun);
04405  
04406   //string temps="myeventsLowEn";
04407   //if (fRun>10000) temps+="MC";
04408   //ofstream& eventsTxtFile=*(this->OpenTxtFile(100,"myevents"));
04409   //ofstream& txtFileLowEn=*(this->OpenTxtFile(100,temps.c_str()));
04410   //ofstream& txtFileNonPC=*(this->OpenTxtFile(100,"myeventsNonPC"));
04411  
04412   //reconstruction object
04413   const MeuReco reco;
04414   const MeuCuts cuts;
04415   const MeuHistos histos;
04416   MeuCounter cnt;
04417 
04418   //variable to turn on all the useful messages if required
04419   Msg::LogLevel_t logLevel=Msg::kDebug;
04420 
04424   
04425   this->InitialiseLoopVariables();  
04426   
04427   //for(Int_t entry=24196;entry<fEntries;entry++){
04428   //for(Int_t entry=20307;entry<fEntries;entry++){//crashed on 32998
04429   //for(Int_t entry=30000;entry<fEntries;entry++){//didn't crash 32998
04430   //for(Int_t entry=0;entry<2000;entry++){
04431   for(Int_t entry=0;entry<fEntries;entry++){
04432     
04433     this->SetLoopVariables(entry);
04434     
04435     MSG("MeuAnalysis",logLevel)
04436       <<"Entry="<<entry<<endl;
04437     
04438     //make a reference instead of a pointer
04439     AtmosEvent& ev=(*fev);
04440     MAXMSG("MeuAnalysis",Msg::kInfo,1)
04441       <<"First entry: run="<<ev.Run
04442       <<", snarl="<<ev.Snarl<<endl;
04443 
04444     //MSG("MeuAnalysis",Msg::kInfo)
04445     //<<"run="<<ev.Run
04446     //<<", snarl="<<ev.Snarl
04447     //<<", entry="<<entry
04448     //<<", unixTime="<<ev.UnixTime
04449     //<<", #scintHits="<<ev.NScintHits
04450     //<<endl;
04451 
04452     //cut out the LI events
04453     if (this->FilterLI(ev)){
04454       //eventsTxtFileLI<<event<<endl;
04455       continue;
04456     }
04457     cnt.notLICounter++;
04458 
04459     //cut out events with multiple tracks
04460     if (ev.NTracks!=1) continue;
04461     cnt.oneTrackCounter++;
04462 
04463     const AtmosTrack* track=dynamic_cast<const AtmosTrack*>
04464       (ev.TrackList->At(0));
04465     const AtmosTrack& trk=(*track);
04466     
04467     //cut out events with less than 20 planes
04468     Bool_t goodLength=trk.AtNuNplanes>20;
04469     if (!goodLength){
04470       continue;
04471     }
04472     cnt.longTrackCounter++;
04473 
04474     //MeuSummary s;
04475     MeuSummary& s=summaryOut.GetMeuSummaryToFill();
04476     s.RFid=3.0*Munits::m;
04477     s.RFidCoil=0.5*Munits::m;
04478     s.DistToEdgeFid=1.0*Munits::m;
04479     s.Event=entry;
04480     cuts.FillSTSumDetails(ev,s);
04481     //a map to hold the info for each plane in the entry
04482     map<Int_t,MeuHitInfo> plInfo;
04483 
04484     cuts.ExtractPlInfo(ev,s,plInfo);
04485     
04486     reco.CalcVtx(s,plInfo);
04487       
04488     MSG("MeuAnalysis",logLevel)<<"CalcFCPC..."<<endl;
04489     reco.CalcFCPC(s,plInfo);
04490 
04491     //check that the XY position is sensible
04492     if (cuts.FilterBadXY(ev,s)){
04493       continue;
04494     }
04495     cnt.goodXYCounter++;
04496     
04497     Bool_t TG=!s.FC && !s.PC;    
04498     if (s.PC) cnt.PCCounter++;
04499     else if (s.FC) cnt.FCCounter++;
04500     else if (TG) cnt.TGCounter++;
04501     if (s.FC && s.PC) MSG("MeuAnalysis",Msg::kError)
04502       <<"Both FC and PC!"<<endl;
04503 
04504     Bool_t awayFromCoil=cuts.AnalyseCoilProximity(ev,s);
04505     if (s.PC && !awayFromCoil) cnt.coilHitCounter++;
04506     if (s.PC && s.SMBoth && awayFromCoil) cnt.bothSMHitCounter++;
04507 
04508     Bool_t goodContainment=s.PC && awayFromCoil && !s.SMBoth;
04509     if (!goodContainment) continue;
04510     MSG("MeuAnalysis",logLevel)<<"Found PC event..."<<endl;
04511     cnt.goodContainmentCounter++;
04512     
04513     Int_t diffAtEnd=-1;
04514     Int_t diffAtVtx=-1;
04515     if (!cuts.CheckTrkViews(s,&diffAtVtx,&diffAtEnd)) {
04516       cnt.badViewsCounter++;
04517       continue;
04518     }
04519 
04520     if (cuts.FilterBadTrackEnd(s,plInfo)){
04521       cnt.badTrackEndCounter++;
04522       continue;
04523     }
04524       
04525     Bool_t goodReco=reco.Reconstruct(s,plInfo);
04526     if (!goodReco){
04527       MSG("MeuAnalysis",logLevel)<<"Bad reco..."<<endl;
04528       cnt.badRecoCounter++;
04529       continue;
04530     }
04531 
04532     if (cuts.FilterBadEndGaps(s,plInfo)){//after reco only
04533       cnt.badEndGapsCounter++;
04534       continue;
04535     }
04536     
04537     MSG("MeuAnalysis",logLevel)<<"Recalibrating..."<<endl;
04538     this->RecalibrateAtNu(ev,s,plInfo);
04539     
04540     MSG("MeuAnalysis",logLevel)<<"CalcWindow..."<<endl;
04541     reco.CalcWindow(s,plInfo);
04542         
04543     //MSG("MeuAnalysis",logLevel)<<"CheckContainment..."<<endl;
04544     //fill the tree once event has been analysed
04545     if (s.WinSigMap<=0) {
04546       cnt.badWindowCounter++;
04547       continue;
04548     }
04549     
04550     if (!reco.CheckWinContainment(s,plInfo)){//nothing for FD
04551       cnt.badFidCounter++;
04552       continue;
04553     }
04554     
04555     //recalculate the containment!
04556     //s.DistToEdgeFid=0.35;//allow a slight shift in position
04557     //reco.CalcFCPC(s,plInfo);
04558     if (!s.PC){
04559       cnt.bad2ndContCounter++;
04560       static Bool_t print=false;
04561       if (print){
04562         MAXMSG("MeuAnalysis",Msg::kInfo,100)
04563           <<endl<<endl<<"Recalculated containment and not PC!"
04564           <<" PC="<<s.PC<<", FC="<<s.FC
04565           <<" entry="<<entry
04566           <<", snarl="<<s.Snarl<<", run="<<s.Run<<endl
04567           <<"Low En="<<s.MCLowEn<<" high="<<s.MCHighEn<<endl
04568           <<endl;
04569         MAXMSG("MeuAnalysis",Msg::kInfo,100)
04570           <<"vtx Pl="<<s.VtxPlane<<", end Pl="<<s.EndPlane
04571           <<", dVtx="<<s.VtxDistToEdge
04572           <<", dEnd="<<s.EndDistToEdge<<endl;
04573         
04574         for (map<Int_t,MeuHitInfo>::const_iterator it=
04575                plInfo.begin();
04576              it!=plInfo.end();++it){
04577           const MeuHitInfo& c=it->second;
04578           MAXMSG("MeuAnalysis",Msg::kInfo,500)
04579             <<"pl="<<c.Plane<<", z="<<c.Z
04580             <<", x="<<c.X<<", y="<<c.Y
04581             <<"  t="<<c.TPos<<", l="<<c.LPos<<endl;
04582         }
04583       }
04584       continue;
04585     }
04586     //txtFileNonPC<<entry<<endl;
04587       
04588     //filter out events with a low traceZ
04589     Float_t trace=999;
04590     Float_t traceZ=999;
04591     reco.CalcTrace(s,plInfo,&trace,&traceZ);
04592     if (traceZ<0.5){
04593       cnt.badTraceZCounter++;
04594       //eventsTxtFileTrace<<event<<endl;
04595       MAXMSG("MeuAnalysis",Msg::kDebug,1000)
04596         <<"trace="<<trace<<", traceZ="<<traceZ
04597         <<", trk.EndTraceZ="<<trk.EndTraceZ<<endl;
04598       continue;
04599     }
04600 
04601     reco.CalcStripDists(s,plInfo);
04602     if (cuts.FilterBadDistEndStrip(s,plInfo)){//nothing for FD
04603       cnt.badDistEndCounter++;
04604       continue;
04605     }
04606     
04608     //MUST BE GOOD!
04610     cnt.goodWindowCounter++;
04611       
04612     cuts.ExtractTruthInfo(ev,s,plInfo);
04613     cuts.FillSTSumRecoDetails(s,plInfo);
04614     histos.FillMeuHistos(s);
04615 
04616     //eventsTxtFile<<entry<<endl;
04617     MAXMSG("MeuAnalysis",Msg::kInfo,10)
04618       <<"Good window: run="
04619       <<ev.Run<<", snarl="<<ev.Snarl<<endl;
04620           
04621     if (s.MCLowEn>0.5*Munits::GeV) {
04622       //txtFileLowEn<<entry<<endl;
04623       MSG("MeuAnalysis",Msg::kInfo)
04624         <<"Low En="<<s.MCLowEn<<" high="<<s.MCHighEn
04625         <<", vtxPl="<<s.VtxPlane<<", endPl="<<s.EndPlane
04626         <<endl<<"entry="<<entry
04627         <<", snarl="<<s.Snarl<<", run="<<s.Run<<endl;
04628     }
04629             
04630     cuts.FillEnVsDist(s,plInfo);
04631     cuts.FillEnVsDist2(s,plInfo);
04632     //cuts.FillTimeHistos(ntp,evt,s);
04633     s.Temperature=this->GetTemperature(s);
04634     summaryOut.SummaryTreeFill(plInfo);
04635   }//end of for                                       
04636     
04640 
04641   MSG("MeuAnalysis",Msg::kInfo)<<"Finished main loop"<<endl;
04642 
04643   MSG("MeuAnalysis",Msg::kInfo) 
04644     <<endl<<"*** Run Summary ***"<<endl;
04645   histos.PrintMeuValues();
04646 
04647   MSG("MeuAnalysis",Msg::kInfo)
04648     <<"Total Entries (snarls)="<<fEntries<<endl;
04649   cnt.PrintNtpSt();
04650 
04651   //finish the summary tree and output it
04652   summaryOut.SummaryTreeFinish();
04653 
04654   MSG("MeuAnalysis",Msg::kInfo) 
04655     <<" ** Finished MakeSummaryTreeWithAtNu method **"<<endl;
04656 }

void MeuAnalysis::MakeSummaryTreeWithNtpSt (  )  [private]

Definition at line 5043 of file MeuAnalysis.cxx.

References fEntries, fRec, fRecBD, InitialiseLoopVariables(), Msg::kInfo, MakeSummaryTreeWithNtpStOneSnarl(), MSG, SetLoopVariables(), and StoreOrFinishSummaryTree().

Referenced by MakeSummaryTree().

05044 {
05045   MSG("MeuAnalysis",Msg::kInfo) 
05046     <<" ** Running MakeSummaryTreeWithNtpSt method... **"<<endl;
05047 
05051   
05052   this->InitialiseLoopVariables();  
05053   
05054   //for(Int_t entry=55000;entry<fEntries;entry++){
05055   for(Int_t entry=0;entry<fEntries;entry++){
05056       
05057     this->SetLoopVariables(entry);
05058 
05059     //run the algorithm to calc meu values for events in snarl
05060     this->MakeSummaryTreeWithNtpStOneSnarl(fRec,fRecBD);
05061 
05062   }//end of for                                       
05063     
05067 
05068   MSG("MeuAnalysis",Msg::kInfo)<<"Finished main loop"<<endl;
05069   
05070   //finish the job
05071   this->StoreOrFinishSummaryTree(NULL,NULL,1);
05072   
05073   MSG("MeuAnalysis",Msg::kInfo) 
05074     <<" ** Finished MakeSummaryTreeWithNtpSt method **"<<endl;
05075 }

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(), Munits::GeV, MeuCounter::goodContainmentCounter, MeuCounter::goodSpillCounter, MeuCounter::goodWindowCounter, MeuCuts::IsAwayFromCoil(), MeuCuts::IsCorrectTrigSrc(), it, Msg::kDebug, Msg::kError, Msg::kInfo, Detector::kNear, Msg::kWarning, 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, NtpStRecord::Print(), print(), MeuCuts::PrintNtpSt(), MeuReco::Reconstruct(), MeuSummary::Run, run(), MeuCuts::SetSpecificCuts(), NtpSREvent::slc, MeuSummary::SMBoth, MeuSummary::Snarl, MeuSummary::SubRun, MeuSummaryWriter::SummaryTreeFill(), MeuSummaryWriter::SummaryTreeSetup(), MeuSummary::Temperature, MeuCounter::TGCounter, MeuHitInfo::TPos, NtpStRecord::trk, NtpSREvent::trk, MeuSummary::VtxDistToEdge, MeuSummary::VtxPlane, MeuSummary::WinSigMap, MeuHitInfo::X, MeuHitInfo::Y, MeuHitInfo::Z, and MeuCounter::zeroEvtsCounter.

Referenced by MeuCalModule::Ana(), and MakeSummaryTreeWithNtpSt().

04662 {
04665 
04666   //get a reference (coding preference)
04667   //don't get a reference to NtpBDLiteRecord because it doesn't exist
04668   //for MC
04669   const NtpStRecord& ntp=*pntp;
04670 
04671   //get candheader
04672   const RecCandHeader& rec=ntp.GetHeader();
04673   Int_t run=rec.GetRun();
04674   Int_t subrun=rec.GetSubRun();
04675   Int_t snarl=rec.GetSnarl();
04676   static Int_t lastRun=-1;
04677   static Int_t lastSubrun=-1;
04678   static Int_t runningTotalOfRuns=0;
04679   static Int_t runningTotalOfSubruns=0;
04680   if (run!=lastRun) {
04681     runningTotalOfRuns++;
04682     MSG("MeuAnalysis",Msg::kInfo)
04683       <<"Found new run="<<run
04684       <<", running total of runs used="<<runningTotalOfRuns<<endl;
04685   }
04686   lastRun=run;
04687   if (subrun!=lastSubrun) {
04688     runningTotalOfSubruns++;
04689     MSG("MeuAnalysis",Msg::kInfo)
04690       <<"Found new subrun="<<subrun
04691       <<", running total of subruns used="<<runningTotalOfSubruns<<endl;
04692   }
04693   lastSubrun=subrun;
04694   MAXMSG("MeuAnalysis",Msg::kInfo,1)
04695     <<"First entry: run="<<run<<", snarl="<<snarl<<endl;
04696   
04697   static const Int_t planeCut=8;
04698   //variable to turn on all the useful messages if required
04699   static const Msg::LogLevel_t logLevel=Msg::kDebug;
04700 
04704   //create static variables to prevent multiple constructions
04705   static Bool_t firstTime=true;
04706 
04707   //reconstruction object
04708   static const MeuReco reco;
04709   static const MeuCuts cuts;
04710   static const MeuHistos histos;
04711   static MeuCounter cnt;
04712   //count all the calls of this loop
04713   cnt.entriesAll++;
04714   
04715   static MeuSummaryWriter* sOut=0;
04716 
04717   static ofstream* peventsTxtFile=0;
04718   static ofstream* ptxtFileLowEn=0;
04719   static ofstream* ptxtFileNonPC=0;
04720 
04721   if (firstTime) {
04722     //set the control variable so code not run again
04723     firstTime=false;
04724     
04725     //create the writer
04726     sOut=new MeuSummaryWriter();
04727     sOut->SummaryTreeSetup(run,subrun);
04728   
04729     //this stores the objects' pointers for access at job end
04730     this->StoreOrFinishSummaryTree(sOut,&cnt,false);
04731 
04732     //printout the first element from the tree
04733     ntp.Print(cout);
04734 
04735     string temps="myeventsLowEn";
04736     if (run>10000) temps+="MC";
04737     peventsTxtFile=this->OpenTxtFile(100,"myevents");
04738     ptxtFileLowEn=this->OpenTxtFile(100,temps.c_str());
04739     ptxtFileNonPC=this->OpenTxtFile(100,"myeventsNonPC");
04740   }
04741 
04742   //get references (coding preference)
04743   MeuSummaryWriter& summaryOut=*sOut;
04744   ofstream& eventsTxtFile=*peventsTxtFile;
04745   ofstream& txtFileLowEn=*ptxtFileLowEn;
04746   ofstream& txtFileNonPC=*ptxtFileNonPC;
04747   
04751 
04752   MSG("MeuAnalysis",logLevel)
04753     <<"Entry="<<cnt.entry<<endl;
04754   
04755   TClonesArray& evtTca=(*ntp.evt);
04756   const Int_t numEvts=evtTca.GetEntriesFast();
04757 
04758   //count the number of snarls with zero events
04759   if (numEvts<=0) cnt.zeroEvtsCounter++;
04760 
04761   //loop over the events
04762   for (Int_t ievt=0;ievt<numEvts;++ievt) {
04763     const NtpSREvent* pevt=
04764       dynamic_cast<NtpSREvent*>(evtTca[ievt]);
04765     const NtpSREvent& evt=(*pevt);
04766 
04767     cnt.evtCounter++;
04768       
04769     //ensure there is only 1 track in the event
04770     if (evt.ntrack!=1) continue;
04771     cnt.oneTrackCounter++;
04772       
04773     //filter out events occuring in a slice with more than 1 event
04774     if (cuts.FilterBadEvtPerSlc(ntp,evt.slc,cnt.entry)) continue;
04775     cnt.oneEvtPerSlcCounter++;
04776 
04777     TClonesArray& trkTca=(*ntp.trk);
04778     const NtpSRTrack* ptrk=
04779       dynamic_cast<NtpSRTrack*>(trkTca[evt.trk[0]]);
04780     const NtpSRTrack& trk=(*ptrk);
04781 
04782     //cut on number of planes
04783     if (sqrt(pow(1.*trk.plane.beg-
04784                  trk.plane.end,2))<planeCut) continue;
04785     cnt.longTrackCounter++;
04786       
04787     //MeuSummary s;
04788     MeuSummary& s=summaryOut.GetMeuSummaryToFill();
04789     s.Event=cnt.entry;
04790     cuts.FillSTSumDetails(ntp,s,ievt,evt.slc);
04791     cuts.SetSpecificCuts(s);
04792     //map to hold the info for each plane in the entry
04793     map<Int_t,MeuHitInfo> plInfo;      
04794       
04795     if (!cuts.IsCorrectTrigSrc(s)){
04796       continue;
04797     }
04798     cnt.correctTrigSrcCounter++;
04799 
04800     //cut out the LI events
04801     if (this->FilterLI(ntp,s)){
04802       //eventsTxtFileLI<<event<<endl;
04803       continue;
04804     }
04805     cnt.notLICounter++;
04806 
04807     if (pntpBD) cuts.GetBDSelectSpillInfo(*pntpBD,s,cnt.entry,
04808                                           cnt.goodSpillCounter);
04809     
04810     //print out info for specific snarl(s)
04811     if (s.Snarl==-1){
04812       MSG("MeuAnalysis",Msg::kInfo)
04813         <<"Entry="<<cnt.entry<<", run="<<run<<", snarl="<<snarl
04814         <<", trk.nstrips="<<trk.nstrip<<endl;
04815       cuts.PrintNtpSt(ntp);
04816     }
04817 
04818     if (cuts.FilterBadTrkTimes(ntp,s,evt)) continue;
04819     cnt.closeTimesCounter++;
04820       
04821     cuts.ExtractPlInfo(ntp,s,plInfo,evt);
04822     if (plInfo.size()==0) {
04823       MAXMSG("MeuAnalysis",Msg::kWarning,20)
04824         <<"No data extracted for entry="<<cnt.entry<<endl;
04825       continue;
04826     }
04827     cnt.notEmptyEventCounter++;
04828       
04829     reco.CalcVtx(s,plInfo);
04830       
04831     //check if track stops in spectrometer
04832     if (s.Detector==Detector::kNear && s.EndPlane>120) continue;
04833     cnt.endPlaneSpectCounter++;
04834 
04835     MSG("MeuAnalysis",logLevel)<<"CalcFCPC..."<<endl;
04836     reco.CalcFCPC(s,plInfo);
04837       
04838     //recalculate the vtx of the tracks
04839     //coming in from the spectrometer
04840     MSG("MeuAnalysis",logLevel)<<"CalcVtxSpecialND..."<<endl;
04841     reco.CalcVtxSpecialND(s,plInfo);
04842 
04843     //cut on number of planes with potential new vertex
04844     if (sqrt(pow(1.*s.VtxPlane-s.EndPlane,2))<planeCut) {
04845       MAXMSG("MeuAnalysis",Msg::kDebug,100)
04846         <<"Cutting on track length, run="<<s.Run
04847         <<", subrun="<<s.SubRun<<", snarl="<<s.Snarl<<endl
04848         <<"s.VtxPlane="<<s.VtxPlane
04849         <<", s.EndPlane="<<s.EndPlane
04850         <<", trk.beg="<<trk.plane.beg
04851         <<", trk.end="<<trk.plane.end<<endl;
04852       //cuts.PrintNtpSt(ntp);
04853       continue;
04854     }
04855     cnt.longTrack2Counter++;
04856       
04857     Bool_t TG=!s.FC && !s.PC;
04858       
04859     if (s.PC) cnt.PCCounter++;
04860     else if (s.FC) cnt.FCCounter++;
04861     else if (TG) cnt.TGCounter++;
04862     if (s.FC && s.PC) MSG("MeuAnalysis",Msg::kError)
04863       <<"Both FC and PC!"<<endl;
04864       
04865     Bool_t awayFromCoil=cuts.IsAwayFromCoil(s,plInfo);
04866     if (s.PC && !awayFromCoil) cnt.coilHitCounter++;
04867     if (s.PC && s.SMBoth && awayFromCoil) cnt.bothSMHitCounter++;
04868       
04869     Bool_t goodContainment=s.PC && !s.SMBoth && awayFromCoil;
04870     if (!goodContainment) continue;
04871     cnt.goodContainmentCounter++;
04872     MSG("MeuAnalysis",logLevel)<<"Found Good PC event..."<<endl;
04873     
04874     Int_t diffAtEnd=-1;
04875     Int_t diffAtVtx=-1;
04876     if (!cuts.CheckTrkViews(s,&diffAtVtx,&diffAtEnd)) {
04877       cnt.badViewsCounter++;
04878       continue;
04879     }
04880 
04881     //check that the hits beyond the end of the track are ok
04882     if (cuts.FilterBadTrackEnd(s,plInfo)){
04883       cnt.badTrackEndCounter++;
04884       continue;
04885     }
04886 
04887     Bool_t goodReco=reco.Reconstruct(s,plInfo);
04888     if (!goodReco){
04889       cnt.badRecoCounter++;
04890       MSG("MeuAnalysis",logLevel)<<"Bad reco..."<<endl;
04891       continue;
04892     }
04893 
04894     //check that there are no big gaps at the end of the track
04895     if (cuts.FilterBadEndGaps(s,plInfo)){//after reco only
04896       cnt.badEndGapsCounter++;
04897       continue;
04898     }
04899 
04900     //recalculate the containment (using MeuReco reconstruction)
04901     s.DistToEdgeFid=0.35;//allow a slight shift in position
04902     reco.CalcFCPC(s,plInfo);
04903           
04904     //re-check containment (using MeuReco reconstruction)
04905     if (!s.PC){
04906       cnt.bad2ndContCounter++;
04907       static Bool_t print=false;
04908       if (print){
04909         MAXMSG("MeuAnalysis",Msg::kInfo,100)
04910           <<endl<<endl<<"Recalculated containment and not PC!"
04911           <<" PC="<<s.PC<<", FC="<<s.FC
04912           <<" entry="<<cnt.entry
04913           <<", snarl="<<s.Snarl<<", run="<<s.Run<<endl
04914           <<"Low En="<<s.MCLowEn<<" high="<<s.MCHighEn<<endl
04915           <<endl;
04916         MAXMSG("MeuAnalysis",Msg::kInfo,100)
04917           <<"vtx Pl="<<s.VtxPlane<<", end Pl="<<s.EndPlane
04918           <<", dVtx="<<s.VtxDistToEdge
04919           <<", dEnd="<<s.EndDistToEdge<<endl;
04920           
04921         for (map<Int_t,MeuHitInfo>::const_iterator it=
04922                plInfo.begin();
04923              it!=plInfo.end();++it){
04924           const MeuHitInfo& c=it->second;
04925           MAXMSG("MeuAnalysis",Msg::kInfo,500)
04926             <<"pl="<<c.Plane<<", z="<<c.Z
04927             <<", x="<<c.X<<", y="<<c.Y
04928             <<"  t="<<c.TPos<<", l="<<c.LPos<<endl;
04929         }          
04930       }
04931       txtFileNonPC<<cnt.entry<<endl;
04932       continue;
04933     }
04934 
04935     MSG("MeuAnalysis",logLevel)<<"Recalibrating..."<<endl;
04936     this->Recalibrate(ntp,s,plInfo,evt);
04937           
04938     MSG("MeuAnalysis",logLevel)<<"CalcWindow..."<<endl;
04939     reco.CalcWindow(s,plInfo);
04940       
04941     MSG("MeuAnalysis",logLevel)<<"CheckContainment..."<<endl;
04942     //check that a track window was found
04943     if (s.WinSigMap<=0) {
04944       cnt.badWindowCounter++;
04945       continue;
04946     }
04947 
04948     //check that the track window is contained within the fid vol
04949     if (!reco.CheckWinContainment(s,plInfo)){//nothing for FD
04950       cnt.badFidCounter++;
04951       continue;
04952     }
04953 
04954     //filter out events with a low traceZ
04955     Float_t trace=999;
04956     Float_t traceZ=999;
04957     reco.CalcTrace(s,plInfo,&trace,&traceZ);//only FD now
04958     if (traceZ<0.5){
04959       cnt.badTraceZCounter++;
04960       //eventsTxtFileTrace<<event<<endl;
04961       MAXMSG("MeuAnalysis",Msg::kDebug,1000)
04962         <<"trace="<<trace<<", traceZ="<<traceZ<<endl;
04963       continue;
04964     }
04965       
04966     //check that all the hits in the track window are >20 cm from
04967     //the end of the strips
04968     reco.CalcStripDists(s,plInfo);
04969     if (cuts.FilterBadDistEndStrip(s,plInfo)){//requires window
04970       cnt.badDistEndCounter++;
04971       continue;
04972     }
04973       
04974     //if beam event check that the track window is away from 
04975     //the possible hadronic shower
04976     if (cuts.FilterBadShwDist(ntp,s,evt)){//only for ND spill data
04977       cnt.badShwDistCounter++;
04978       continue;
04979     }
04980       
04982     //MUST BE GOOD!
04984     cnt.goodWindowCounter++;
04985             
04986     cuts.ExtractTruthInfo(ntp,s,plInfo);
04987     cuts.FillSTSumRecoDetails(s,plInfo);
04988     histos.FillMeuHistos(s);
04989     if (pntpBD) cuts.FillBeamMonDetails(*pntpBD,s);
04990     
04991     eventsTxtFile<<cnt.entry<<endl;
04992     MAXMSG("MeuAnalysis",Msg::kInfo,10)
04993       <<"Good track window: run="<<run<<", snarl="<<snarl
04994       <<", evt="<<ievt+1<<"/"<<numEvts<<endl;
04995       
04996     if (s.MCLowEn>0.5*Munits::GeV) {
04997       txtFileLowEn<<cnt.entry<<endl;
04998       MSG("MeuAnalysis",Msg::kInfo)
04999         <<"Low En="<<s.MCLowEn<<" high="<<s.MCHighEn
05000         <<", vtxPl="<<s.VtxPlane<<", endPl="<<s.EndPlane
05001         <<endl<<"entry="<<cnt.entry
05002         <<", snarl="<<s.Snarl<<", run="<<s.Run<<endl;
05003     }
05004       
05005     cuts.FillEnVsDist(s,plInfo);
05006     cuts.FillEnVsDist2(s,plInfo);
05007     cuts.FillTimeHistos(ntp,evt,s);//a loop over lots of strips...
05008     s.Temperature=this->GetTemperature(s);
05009     summaryOut.SummaryTreeFill(plInfo);
05010   }//end of loop over events
05011 
05012   //increment ntuple entry (snarl) counter
05013   cnt.entry++;
05014 }

void MeuAnalysis::N_1Plots (  ) 

Definition at line 1114 of file MeuAnalysis.cxx.

References 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(), fEntries, MeuCuts::FillSTSumDetails(), MeuCuts::FilterBadDistEndStrip(), MeuCuts::FilterBadEndGaps(), MeuCuts::FilterBadTrackEnd(), MeuCuts::FilterLowMatTraversed(), fOutFile, InitialiseLoopVariables(), Msg::kInfo, MAXMSG, MSG, OpenFile(), MeuSummary::PC, NtpSRVertex::plane, Recalibrate(), MeuReco::Reconstruct(), MeuSummary::Reset(), MeuSummary::RFid, SetLoopVariables(), MeuSummary::TotalMatTraversed, NtpStRecord::trk, NtpSRTrack::vtx, MeuSummary::VtxPlane, MeuSummary::WinAvPLCor, MeuSummary::WinSigMap, MeuSummary::WinStopSidePl, and MeuSummary::WinVtxSidePl.

01115 {
01116   MSG("MeuAnalysis",Msg::kInfo) 
01117     <<" ** Running N_1Plots method... **"<<endl;
01118 
01119   //open the output file for the histograms
01120   fOutFile=this->OpenFile(100,"N_1Plots");
01121 
01122   TH1F* hTrkPl=new TH1F("hTrkPl","hTrkPl",201,-1,200);
01123   hTrkPl->SetTitle("hTrkPl");
01124   hTrkPl->GetXaxis()->SetTitle("hTrkPl");
01125   hTrkPl->GetXaxis()->CenterTitle();
01126   hTrkPl->SetFillColor(0);
01127   hTrkPl->SetLineWidth(3);
01128   //hTrkPl->SetBit(TH1::kCanRebin);
01129 
01130   TH1F* hTrkPlN_1=new TH1F("hTrkPlN_1","hTrkPlN_1",201,-1,200);
01131   hTrkPlN_1->SetTitle("hTrkPlN_1");
01132   hTrkPlN_1->GetXaxis()->SetTitle("hTrkPlN_1");
01133   hTrkPlN_1->GetXaxis()->CenterTitle();
01134   hTrkPlN_1->SetLineColor(2);//red for N-1 cuts
01135   hTrkPlN_1->SetFillColor(0);
01136   hTrkPlN_1->SetLineWidth(3);
01137   //hTrkPlN_1->SetBit(TH1::kCanRebin);
01138 
01139   TH1F* hTrkPlN=new TH1F("hTrkPlN","hTrkPlN",201,-1,200);
01140   hTrkPlN->SetTitle("hTrkPlN");
01141   hTrkPlN->GetXaxis()->SetTitle("hTrkPlN");
01142   hTrkPlN->GetXaxis()->CenterTitle();
01143   hTrkPlN->SetLineColor(4);//blue for N cuts
01144   hTrkPlN->SetFillColor(0);
01145   hTrkPlN->SetLineWidth(3);
01146   //hTrkPlN->SetBit(TH1::kCanRebin);
01147 
01148   TH1F* hDistEdge=new TH1F("hDistEdge","hDistEdge",201,-3,3);
01149   hDistEdge->SetTitle("hDistEdge");
01150   hDistEdge->GetXaxis()->SetTitle("hDistEdge");
01151   hDistEdge->GetXaxis()->CenterTitle();
01152   hDistEdge->SetFillColor(0);
01153   hDistEdge->SetLineWidth(3);
01154   //hDistEdge->SetBit(TH1::kCanRebin);
01155 
01156   TH1F* hDistEdgeN_1=new TH1F("hDistEdgeN_1","hDistEdgeN_1",201,-3,3);
01157   hDistEdgeN_1->SetTitle("hDistEdgeN_1");
01158   hDistEdgeN_1->GetXaxis()->SetTitle("hDistEdgeN_1");
01159   hDistEdgeN_1->GetXaxis()->CenterTitle();
01160   hDistEdgeN_1->SetLineColor(2);//red for N-1 cuts
01161   hDistEdgeN_1->SetFillColor(0);
01162   hDistEdgeN_1->SetLineWidth(3);
01163   //hDistEdgeN_1->SetBit(TH1::kCanRebin);
01164 
01165   TH1F* hDistEdgeN=new TH1F("hDistEdgeN","hDistEdgeN",201,-3,3);
01166   hDistEdgeN->SetTitle("hDistEdgeN");
01167   hDistEdgeN->GetXaxis()->SetTitle("hDistEdgeN");
01168   hDistEdgeN->GetXaxis()->CenterTitle();
01169   hDistEdgeN->SetLineColor(4);//blue for N cuts
01170   hDistEdgeN->SetFillColor(0);
01171   hDistEdgeN->SetLineWidth(3);
01172   //hDistEdgeN->SetBit(TH1::kCanRebin);
01173 
01174   TH1F* hViewsVtx=new TH1F("hViewsVtx","hViewsVtx",101,-1,100);
01175   hViewsVtx->SetTitle("hViewsVtx");
01176   hViewsVtx->GetXaxis()->SetTitle("hViewsVtx");
01177   hViewsVtx->GetXaxis()->CenterTitle();
01178   hViewsVtx->SetFillColor(0);
01179   hViewsVtx->SetLineWidth(3);
01180   //hViewsVtx->SetBit(TH1::kCanRebin);
01181 
01182   TH1F* hViewsVtxN_1=new TH1F("hViewsVtxN_1","hViewsVtxN_1",101,-1,100);
01183   hViewsVtxN_1->SetTitle("hViewsVtxN_1");
01184   hViewsVtxN_1->GetXaxis()->SetTitle("hViewsVtxN_1");
01185   hViewsVtxN_1->GetXaxis()->CenterTitle();
01186   hViewsVtxN_1->SetLineColor(2);//red for N-1 cuts
01187   hViewsVtxN_1->SetFillColor(0);
01188   hViewsVtxN_1->SetLineWidth(3);
01189   //hViewsVtxN_1->SetBit(TH1::kCanRebin);
01190 
01191   TH1F* hViewsVtxN=new TH1F("hViewsVtxN","hViewsVtxN",101,-1,100);
01192   hViewsVtxN->SetTitle("hViewsVtxN");
01193   hViewsVtxN->GetXaxis()->SetTitle("hViewsVtxN");
01194   hViewsVtxN->GetXaxis()->CenterTitle();
01195   hViewsVtxN->SetLineColor(4);//blue for N cuts
01196   hViewsVtxN->SetFillColor(0);
01197   hViewsVtxN->SetLineWidth(3);
01198   //hViewsVtxN->SetBit(TH1::kCanRebin);
01199 
01200   TH1F* hViewsEnd=new TH1F("hViewsEnd","hViewsEnd",101,-1,100);
01201   hViewsEnd->SetTitle("hViewsEnd");
01202   hViewsEnd->GetXaxis()->SetTitle("hViewsEnd");
01203   hViewsEnd->GetXaxis()->CenterTitle();
01204   hViewsEnd->SetFillColor(0);
01205   hViewsEnd->SetLineWidth(3);
01206   //hViewsEnd->SetBit(TH1::kCanRebin);
01207 
01208   TH1F* hViewsEndN_1=new TH1F("hViewsEndN_1","hViewsEndN_1",101,-1,100);
01209   hViewsEndN_1->SetTitle("hViewsEndN_1");
01210   hViewsEndN_1->GetXaxis()->SetTitle("hViewsEndN_1");
01211   hViewsEndN_1->GetXaxis()->CenterTitle();
01212   hViewsEndN_1->SetLineColor(2);//red for N-1 cuts
01213   hViewsEndN_1->SetFillColor(0);
01214   hViewsEndN_1->SetLineWidth(3);
01215   //hViewsEndN_1->SetBit(TH1::kCanRebin);
01216 
01217   TH1F* hViewsEndN=new TH1F("hViewsEndN","hViewsEndN",101,-1,100);
01218   hViewsEndN->SetTitle("hViewsEndN");
01219   hViewsEndN->GetXaxis()->SetTitle("hViewsEndN");
01220   hViewsEndN->GetXaxis()->CenterTitle();
01221   hViewsEndN->SetLineColor(4);//blue for N cuts
01222   hViewsEndN->SetFillColor(0);
01223   hViewsEndN->SetLineWidth(3);
01224   //hViewsEndN->SetBit(TH1::kCanRebin);
01225 
01226   TH1F* hDistEndSt=new TH1F("hDistEndSt","hDistEndSt",200,-0.1,2);
01227   hDistEndSt->SetTitle("hDistEndSt");
01228   hDistEndSt->GetXaxis()->SetTitle("hDistEndSt");
01229   hDistEndSt->GetXaxis()->CenterTitle();
01230   hDistEndSt->SetFillColor(0);
01231   hDistEndSt->SetLineWidth(3);
01232   //hDistEndSt->SetBit(TH1::kCanRebin);
01233 
01234   TH1F* hDistEndStN_1=new TH1F("hDistEndStN_1","hDistEndStN_1",
01235                                200,-0.1,2);
01236   hDistEndStN_1->SetTitle("hDistEndStN_1");
01237   hDistEndStN_1->GetXaxis()->SetTitle("hDistEndStN_1");
01238   hDistEndStN_1->GetXaxis()->CenterTitle();
01239   hDistEndStN_1->SetLineColor(2);//red for N-1 cuts
01240   hDistEndStN_1->SetFillColor(0);
01241   hDistEndStN_1->SetLineWidth(3);
01242   //hDistEndStN_1->SetBit(TH1::kCanRebin);
01243 
01244   TH1F* hDistEndStN=new TH1F("hDistEndStN","hDistEndStN",200,-0.1,2);
01245   hDistEndStN->SetTitle("hDistEndStN");
01246   hDistEndStN->GetXaxis()->SetTitle("hDistEndStN");
01247   hDistEndStN->GetXaxis()->CenterTitle();
01248   hDistEndStN->SetLineColor(4);//blue for N cuts
01249   hDistEndStN->SetFillColor(0);
01250   hDistEndStN->SetLineWidth(3);
01251   //hDistEndStN->SetBit(TH1::kCanRebin);
01252 
01253   TH1F* hSigBeyondEnd=new TH1F("hSigBeyondEnd","hSigBeyondEnd",1000,-1,50000);
01254   hSigBeyondEnd->SetTitle("hSigBeyondEnd");
01255   hSigBeyondEnd->GetXaxis()->SetTitle("hSigBeyondEnd");
01256   hSigBeyondEnd->GetXaxis()->CenterTitle();
01257   hSigBeyondEnd->SetFillColor(0);
01258   hSigBeyondEnd->SetLineWidth(3);
01259   //hSigBeyondEnd->SetBit(TH1::kCanRebin);
01260   
01261   TH1F* hSigBeyondEndN_1=new TH1F("hSigBeyondEndN_1","hSigBeyondEndN_1",
01262                                   1000,-1,50000);
01263   hSigBeyondEndN_1->SetTitle("hSigBeyondEndN_1");
01264   hSigBeyondEndN_1->GetXaxis()->SetTitle("hSigBeyondEndN_1");
01265   hSigBeyondEndN_1->GetXaxis()->CenterTitle();
01266   hSigBeyondEndN_1->SetLineColor(2);//red for N-1 cuts
01267   hSigBeyondEndN_1->SetFillColor(0);
01268   hSigBeyondEndN_1->SetLineWidth(3);
01269   //hSigBeyondEndN_1->SetBit(TH1::kCanRebin);
01270 
01271   TH1F* hSigBeyondEndN=new TH1F("hSigBeyondEndN","hSigBeyondEndN",
01272                                 1000,-1,50000);
01273   hSigBeyondEndN->SetTitle("hSigBeyondEndN");
01274   hSigBeyondEndN->GetXaxis()->SetTitle("hSigBeyondEndN");
01275   hSigBeyondEndN->GetXaxis()->CenterTitle();
01276   hSigBeyondEndN->SetLineColor(4);//blue for N cuts
01277   hSigBeyondEndN->SetFillColor(0);
01278   hSigBeyondEndN->SetLineWidth(3);
01279   //hSigBeyondEndN->SetBit(TH1::kCanRebin);
01280 
01281 
01282   TH1F* hSigTrkEndGaps=new TH1F("hSigTrkEndGaps","hSigTrkEndGaps",1000,-1,50000);
01283   hSigTrkEndGaps->SetTitle("hSigTrkEndGaps");
01284   hSigTrkEndGaps->GetXaxis()->SetTitle("hSigTrkEndGaps");
01285   hSigTrkEndGaps->GetXaxis()->CenterTitle();
01286   hSigTrkEndGaps->SetFillColor(0);
01287   hSigTrkEndGaps->SetLineWidth(3);
01288   //hSigTrkEndGaps->SetBit(TH1::kCanRebin);
01289   
01290   TH1F* hSigTrkEndGapsN_1=new TH1F("hSigTrkEndGapsN_1","hSigTrkEndGapsN_1",
01291                                   1000,-1,50000);
01292   hSigTrkEndGapsN_1->SetTitle("hSigTrkEndGapsN_1");
01293   hSigTrkEndGapsN_1->GetXaxis()->SetTitle("hSigTrkEndGapsN_1");
01294   hSigTrkEndGapsN_1->GetXaxis()->CenterTitle();
01295   hSigTrkEndGapsN_1->SetLineColor(2);//red for N-1 cuts
01296   hSigTrkEndGapsN_1->SetFillColor(0);
01297   hSigTrkEndGapsN_1->SetLineWidth(3);
01298   //hSigTrkEndGapsN_1->SetBit(TH1::kCanRebin);
01299 
01300   TH1F* hSigTrkEndGapsN=new TH1F("hSigTrkEndGapsN","hSigTrkEndGapsN",
01301                                 1000,-1,50000);
01302   hSigTrkEndGapsN->SetTitle("hSigTrkEndGapsN");
01303   hSigTrkEndGapsN->GetXaxis()->SetTitle("hSigTrkEndGapsN");
01304   hSigTrkEndGapsN->GetXaxis()->CenterTitle();
01305   hSigTrkEndGapsN->SetLineColor(4);//blue for N cuts
01306   hSigTrkEndGapsN->SetFillColor(0);
01307   hSigTrkEndGapsN->SetLineWidth(3);
01308   //hSigTrkEndGapsN->SetBit(TH1::kCanRebin);
01309 
01310   TH1F* hMatTrav=new TH1F("hMatTrav","hMatTrav",400,0,10);
01311   hMatTrav->SetTitle("hMatTrav");
01312   hMatTrav->GetXaxis()->SetTitle("hMatTrav");
01313   hMatTrav->GetXaxis()->CenterTitle();
01314   hMatTrav->SetFillColor(0);
01315   hMatTrav->SetLineWidth(3);
01316   //hMatTrav->SetBit(TH1::kCanRebin);
01317 
01318   TH1F* hMatTravN_1=new TH1F("hMatTravN_1","hMatTravN_1",400,0,10);
01319   hMatTravN_1->SetTitle("hMatTravN_1");
01320   hMatTravN_1->GetXaxis()->SetTitle("hMatTravN_1");
01321   hMatTravN_1->GetXaxis()->CenterTitle();
01322   hMatTravN_1->SetLineColor(2);//red for N-1 cuts
01323   hMatTravN_1->SetFillColor(0);
01324   hMatTravN_1->SetLineWidth(3);
01325   //hMatTravN_1->SetBit(TH1::kCanRebin);
01326 
01327   TH1F* hMatTravN=new TH1F("hMatTravN","hMatTravN",400,0,10);
01328   hMatTravN->SetTitle("hMatTravN");
01329   hMatTravN->GetXaxis()->SetTitle("hMatTravN");
01330   hMatTravN->GetXaxis()->CenterTitle();
01331   hMatTravN->SetLineColor(4);//blue for N cuts
01332   hMatTravN->SetFillColor(0);
01333   hMatTravN->SetLineWidth(3);
01334   //hMatTravN->SetBit(TH1::kCanRebin);
01335         
01336   //reconstruction object
01337   MeuReco reco;
01338   MeuCuts cuts;
01339   MeuSummary s;
01340 
01344   
01345   this->InitialiseLoopVariables();  
01346   
01347   //for(Int_t entry=0;entry<10000;entry++){
01348   for(Int_t entry=0;entry<fEntries;entry++){
01349       
01350     this->SetLoopVariables(entry);
01351       
01352     NtpStRecord& ntp=(*fRec);//Make a reference instead of a pointer
01353 
01354     //ensure there is only 1 track in the entry
01355     TClonesArray& trkTca=(*ntp.trk);
01356     Int_t numTrks=trkTca.GetEntries();
01357     if (numTrks!=1) continue;
01358     const NtpSRTrack* ptrk=dynamic_cast<NtpSRTrack*>(trkTca[0]);
01359     const NtpSRTrack& trk=(*ptrk);
01360    
01361     //cut on number of planes
01362     Int_t numPlanes=abs(trk.vtx.plane-trk.end.plane);
01363     if (numPlanes<8) continue;
01364 
01365     s.Reset();
01366     s.RFid=1.0;
01367     s.DistToEdgeFid=0.4;
01368     s.Event=entry;
01369     cuts.FillSTSumDetails(ntp,s,-1,-1);
01370     map<Int_t,MeuHitInfo> plInfo;
01371 
01372     //check if the track is in the spectrometer
01373     if (trk.vtx.plane>120 && trk.end.plane>120) continue;
01374 
01375     TClonesArray& evtTca=(*ntp.evt);
01376     const NtpSREvent* pevt=dynamic_cast<NtpSREvent*>(evtTca[0]);
01377     cuts.ExtractPlInfo(ntp,s,plInfo,*pevt);
01378     
01379     reco.CalcVtx(s,plInfo);
01380 
01381     //check if track stops in spectrometer
01382     if (s.EndPlane>120) continue;
01383 
01384     //calculate the containment (before doing the vtxSpecial)
01385     reco.CalcFCPC(s,plInfo);
01386 
01387     //recalculate the vtx of the tracks coming in from the spectrometer
01388     reco.CalcVtxSpecialND(s,plInfo);
01389 
01390     //cut on number of planes with potential new vertex
01391     
01392     numPlanes=abs(s.VtxPlane-s.EndPlane);
01393     if (numPlanes<8) continue;
01394 
01395     reco.CalcStripDists(s,plInfo);
01396     
01397     //calc distance to detector edge
01398     Float_t dist=reco.CalcSmallestDeepDistToEdge(plInfo[s.EndPlane]);
01399 
01400     Int_t diffAtEnd=-1;
01401     Int_t diffAtVtx=-1;
01402     cuts.CheckTrkViews(s,&diffAtVtx,&diffAtEnd);
01403 
01404     Float_t minDist=-1;
01405     Float_t sigBeyondEnd=-1;
01406     Float_t adcInGap=-1;
01407     
01408     Bool_t goodContainment=false;
01409     Bool_t goodReco=false;
01410     if (s.PC) {
01411       goodReco=reco.Reconstruct(s,plInfo);
01412       if (goodReco) {
01413 
01414         TClonesArray& evtTca=(*ntp.evt);
01415         const NtpSREvent* pevt=dynamic_cast<NtpSREvent*>(evtTca[0]);
01416         this->Recalibrate(ntp,s,plInfo,*pevt);
01417         reco.CalcWindow(s,plInfo);
01418 
01419         //now check that window is deeply contained
01420         if (s.WinSigMap>0) {
01421           goodContainment=reco.CheckWinContainment(s,plInfo);
01422         }
01423 
01424         //now do some of the calculations that require good PC
01425         if (goodContainment){
01426           cuts.FilterBadDistEndStrip(s,plInfo,&minDist);
01427           cuts.FilterBadTrackEnd(s,plInfo,&sigBeyondEnd);
01428           cuts.FilterBadEndGaps(s,plInfo,&adcInGap);
01429         }
01430       }
01431     }
01432 
01433 
01434     //if (s.WinSigMap<=0) { //matTrav
01435     
01436     
01440     //start of "no cuts"
01441     hTrkPl->Fill(numPlanes);
01442     hDistEdge->Fill(dist);
01443     hViewsVtx->Fill(diffAtVtx);
01444     hViewsEnd->Fill(diffAtEnd);
01445     if (goodContainment){ //only have a value in this case
01446       hDistEndSt->Fill(minDist);
01447       hSigBeyondEnd->Fill(sigBeyondEnd);
01448       hSigTrkEndGaps->Fill(adcInGap);
01449     }
01450     //only have a value in this case (requires calcWindow)
01451     if (goodReco) hMatTrav->Fill(s.TotalMatTraversed);
01452     
01453     //cuts
01454     Bool_t goodDist=dist>0.4;
01455     Bool_t goodTrkPl=numPlanes>10;
01456     Bool_t goodViewsVtx=diffAtVtx<=9;//for ND full inst. planes
01457     Bool_t goodViewsEnd=diffAtEnd<=5;
01458     Bool_t goodDistEndSt=minDist>0.2;
01459     Bool_t goodSigBeyondEnd=sigBeyondEnd<500 && sigBeyondEnd>=0;
01460     Bool_t goodSigTrkEndGaps=adcInGap<500 && adcInGap>=0;
01461     Bool_t goodMatTrav=!cuts.FilterLowMatTraversed(s);
01462   
01463     Bool_t goodN=
01464       goodContainment && 
01465       goodDist && 
01466       goodTrkPl && 
01467       goodViewsVtx &&
01468       goodViewsEnd &&
01469       goodDistEndSt &&
01470       goodSigBeyondEnd &&
01471       goodSigTrkEndGaps &&
01472       goodMatTrav &&
01473       true;//for coding convenience
01474     
01478     //start of N-1 cuts plots
01479     if (goodContainment && 
01480         goodDist &&
01481         goodViewsVtx &&
01482         goodViewsEnd &&
01483         goodDistEndSt &&
01484         goodSigBeyondEnd &&
01485         goodSigTrkEndGaps &&
01486         goodMatTrav &&
01487         true) {
01488       hTrkPlN_1->Fill(numPlanes);
01489 
01490       if (numPlanes<10){
01491         MAXMSG("MeuAnalysis",Msg::kInfo,10)
01492           <<"vtxPl="<<s.VtxPlane
01493           <<", endPl="<<s.EndPlane
01494           <<", stopPl="<<s.WinVtxSidePl
01495           <<", startPl="<<s.WinStopSidePl
01496           <<", meu="<<s.WinSigMap
01497           <<", PLC="<<s.WinAvPLCor
01498           <<endl;
01499       }
01500     }
01501     if (goodContainment && 
01502         goodTrkPl &&
01503         goodViewsVtx &&
01504         goodViewsEnd &&
01505         goodDistEndSt &&
01506         goodSigBeyondEnd &&
01507         goodSigTrkEndGaps &&
01508         goodMatTrav &&
01509         true) {
01510       hDistEdgeN_1->Fill(dist);
01511     }
01512     if (goodContainment && 
01513         goodDist &&
01514         goodTrkPl &&
01515         goodViewsEnd &&
01516         goodDistEndSt &&
01517         goodSigBeyondEnd &&
01518         goodSigTrkEndGaps &&
01519         goodMatTrav &&
01520         true) {
01521       hViewsVtxN_1->Fill(diffAtVtx);
01522     }
01523     if (goodContainment && 
01524         goodDist &&
01525         goodTrkPl &&
01526         goodViewsVtx &&
01527         goodDistEndSt &&
01528         goodSigBeyondEnd &&
01529         goodSigTrkEndGaps &&
01530         goodMatTrav &&
01531         true) {
01532       hViewsEndN_1->Fill(diffAtEnd);
01533     }
01534     if (goodContainment && 
01535         goodDist &&
01536         goodTrkPl &&
01537         goodViewsVtx &&
01538         goodViewsEnd &&
01539         goodSigBeyondEnd &&
01540         goodSigTrkEndGaps &&
01541         goodMatTrav &&
01542         true) {
01543       hDistEndStN_1->Fill(minDist);
01544     }
01545     if (goodContainment && 
01546         goodDist &&
01547         goodTrkPl &&
01548         goodViewsVtx &&
01549         goodViewsEnd &&
01550         goodDistEndSt &&
01551         goodSigTrkEndGaps &&
01552         goodMatTrav &&
01553         true) {
01554       hSigBeyondEndN_1->Fill(sigBeyondEnd);
01555     }
01556     if (goodContainment && 
01557         goodDist &&
01558         goodTrkPl &&
01559         goodViewsVtx &&
01560         goodViewsEnd &&
01561         goodDistEndSt &&
01562         goodSigBeyondEnd &&
01563         goodMatTrav &&
01564         true) {
01565       hSigTrkEndGapsN_1->Fill(adcInGap);
01566     }
01567     if (goodContainment && 
01568         goodDist &&
01569         goodTrkPl &&
01570         goodViewsVtx &&
01571         goodViewsEnd &&
01572         goodDistEndSt &&
01573         goodSigBeyondEnd &&
01574         goodSigTrkEndGaps &&
01575         true) {
01576       hMatTravN_1->Fill(s.TotalMatTraversed);
01577     }    
01578 
01582     //start of N cuts plots
01583     if (goodN) {
01584       hTrkPlN->Fill(numPlanes);
01585       hDistEdgeN->Fill(dist);
01586       hViewsVtxN->Fill(diffAtVtx);
01587       hViewsEndN->Fill(diffAtEnd);
01588       hDistEndStN->Fill(minDist);
01589       hSigBeyondEndN->Fill(sigBeyondEnd);
01590       hSigTrkEndGapsN->Fill(adcInGap);
01591       hMatTravN->Fill(s.TotalMatTraversed);
01592     }
01593   }//end of for                                       
01594   
01598 
01599   MSG("MeuAnalysis",Msg::kInfo)<<"Finished main loop"<<endl;
01600 
01601   MSG("MeuAnalysis",Msg::kInfo) 
01602     <<" ** Finished N_1Plots method **"<<endl;
01603 }

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

Definition at line 338 of file MeuAnalysis.cxx.

References Msg::kInfo, Msg::kWarning, and MSG.

Referenced by MakeChain().

00340 {
00341   if (!file) {
00342     MSG("MeuAnalysis",Msg::kWarning)
00343       <<"File does not exist so can't test if object exists!"<<endl;
00344     return false;
00345   }
00346 
00347   TList* listOfKeys=file->GetListOfKeys();
00348 
00349   static TFile* lastFile=0;
00350   if (lastFile!=file){
00351     MSG("MeuAnalysis",Msg::kInfo)
00352       <<"File contains these keys:"<<endl;
00353     listOfKeys->Print();
00354   }
00355   lastFile=file;
00356 
00357   TObject* ob=listOfKeys->FindObject(objectName.c_str());
00358   if (ob) {
00359     MSG("MeuAnalysis",Msg::kInfo)
00360       <<"Requested object with name="<<objectName
00361       <<" exists in file:"<<endl;
00362     ob->Print();
00363     return true;
00364   }
00365   else {
00366     MSG("MeuAnalysis",Msg::kInfo)
00367       <<"Requested object with name="<<objectName
00368       <<" does NOT exist in file"<<endl;
00369     return false;
00370   }
00371 }

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

Referenced by BasicPlots(), N_1Plots(), and SpillPlots().

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

Definition at line 712 of file MeuAnalysis.cxx.

References Form(), Msg::kInfo, and MSG.

Referenced by BasicReco(), SnarlList(), and SpillPlots().

00714 {
00715   //create the tfile pointer to be returned
00716   ofstream* outputFile=0;
00717   
00718   //convert variables to string
00719   string sRunNumber=Form("%d",runNumber);
00720   string sDetector="F";
00721   string sPrefix="";//default
00722   string sAnaDir=".";
00723   if (prefix!="") sPrefix+=prefix;
00724   string sBase=sAnaDir+"/"+sPrefix+sDetector+sRunNumber;
00725   string sFileName=sBase+".txt";
00726 
00727   outputFile=new ofstream(sFileName.c_str());
00728   
00729   if (outputFile){
00730     MSG("MeuAnalysis",Msg::kInfo) 
00731       <<"Output txt file opened: "<<sFileName<<endl;
00732   }
00733   else{
00734     MSG("MeuAnalysis",Msg::kInfo)
00735       <<"Txt file NOT opened: "<<sFileName<<endl;
00736   }
00737 
00738   return outputFile;
00739 }

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, Detector::kFar, Detector::kNear, Msg::kVerbose, Msg::kWarning, and MSG.

Referenced by Recalibrate().

02049 {
02050   Bool_t withinTrk=false;
02051   
02052   if (s.Detector==Detector::kNear){
02053     if (plane<0 || plane>120) {
02054       MSG("MeuAnalysis",Msg::kVerbose)
02055         <<"MeuAnalysis::PlaneIsWithinTrack."
02056         <<" Plane is out of range="<<plane
02057         <<", det="<<s.Detector<<endl;
02058       return withinTrk;
02059     }
02060   }
02061   else if (s.Detector==Detector::kFar){
02062     if (plane<0 || plane>485) {
02063       MSG("MeuAnalysis",Msg::kVerbose)
02064         <<"MeuAnalysis::PlaneIsWithinTrack."
02065         <<" Plane is out of range="<<plane
02066         <<", det="<<s.Detector<<endl;
02067       return withinTrk;
02068     }
02069   }
02070   else {
02071     MSG("MeuAnalysis",Msg::kWarning)
02072       <<"MeuAnalysis::PlaneIsWithinTrack."
02073       <<" Detector not supported"<<plane
02074       <<endl; 
02075   }
02076 
02077 
02078   //count the end planes as being "within" the track
02079   if ((plane>=vtxPl && plane<=endPl) ||
02080       (plane<=vtxPl && plane>=endPl)) withinTrk=true;
02081   
02082   return withinTrk;
02083 }

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, fRecalibratePE, fRecalibrateSigCor, fRecalibrateSigLin, Calibrator::GetAttenCorrectedTpos(), Calibrator::GetDriftCorrected(), Calibrator::GetLinearized(), Calibrator::GetPhotoElectrons(), Calibrator::GetStripToStripCorrected(), Calibrator::Instance(), it, Msg::kDebug, StripEnd::kEast, Msg::kInfo, SimFlag::kMC, Detector::kNear, StripEnd::kWest, MeuHitInfo::LPos, MAXMSG, MeuHitInfo::MCLPos, MeuHitInfo::MCSigMap, MeuSummary::MedianTime, Munits::ns, NtpSRSlice::nstrip, MeuHitInfo::Pe, NtpSRPulseHeight::pe, MeuHitInfo::Pe1, MeuHitInfo::Pe2, NtpSRStrip::ph0, NtpSRStrip::ph1, MeuHitInfo::Plane, NtpSRStrip::plane, PlaneIsWithinTrack(), print(), NtpSRPulseHeight::raw, Calibrator::ReInitialise(), 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, NtpStRecord::stp, NtpSRSlice::stp, NtpSRStrip::strip, MeuHitInfo::Strip, NtpSRStrip::time0, NtpSRStrip::time1, MeuSummary::TimeSec, and MeuSummary::VtxPlane.

Referenced by BasicReco(), N_1Plots(), and SpillPlots().

02129 {
02132   
02133   //variable to turn on all the useful messages if required
02134   Msg::LogLevel_t logLevel=Msg::kDebug;
02135   //Msg::LogLevel_t logLevelSts=Msg::kDebug;
02136   //Msg::LogLevel_t logLevelLin=Msg::kInfo;
02137 
02138   //make plots of the calibration constants for the recalibration
02139   //this is also done elsewhere for ntuple values of the cal constants
02140   //if the recalibration is not done then they are the same as sntp
02141   static TH1F* hAdcPerSigLinReCal1=0;
02142   static TH1F* hSigLinPerSigCorReCal1=0;
02143   static TH1F* hSigCorPerSigMapReCal1=0;
02144   static TH1F* hSigMapPerMipReCal1=0;
02145   static TH1F* hMipPerGeVReCal1=0;
02146   static TH1F* hAdcPerPeReCal1=0;
02147   static TH1F* hAdcPerSigDrfReCal1=0;
02148   static TH1F* hAdcPerSigLinOnlyReCal1=0;
02149 
02150   static TH1F* hAdcPerSigLinReCal2=0;
02151   static TH1F* hSigLinPerSigCorReCal2=0;
02152   static TH1F* hSigCorPerSigMapReCal2=0;
02153   static TH1F* hSigMapPerMipReCal2=0;
02154   static TH1F* hMipPerGeVReCal2=0;
02155   static TH1F* hAdcPerPeReCal2=0;
02156   static TH1F* hAdcPerSigDrfReCal2=0;
02157   static TH1F* hAdcPerSigLinOnlyReCal2=0;
02158 
02159   if (hAdcPerSigLinReCal1==0) {
02160     //stripend 1
02161     hAdcPerSigLinReCal1=new TH1F("hAdcPerSigLinReCal1",
02162                                  "hAdcPerSigLinReCal1",
02163                                  1100,-1,10);
02164     hAdcPerSigLinReCal1->SetFillColor(0);
02165     hAdcPerSigLinReCal1->SetTitle("AdcPerSigLin");
02166     hAdcPerSigLinReCal1->GetXaxis()->SetTitle("AdcPerSigLin");
02167     hAdcPerSigLinReCal1->GetXaxis()->CenterTitle();
02168     //hAdcPerSigLinReCal1->SetBit(TH1::kCanRebin);
02169     
02170     hSigLinPerSigCorReCal1=new TH1F("hSigLinPerSigCorReCal1",
02171                                  "hSigLinPerSigCorReCal1",
02172                                  1100,-1,10);
02173     hSigLinPerSigCorReCal1->SetFillColor(0);
02174     hSigLinPerSigCorReCal1->SetTitle("SigLinPerSigCor");
02175     hSigLinPerSigCorReCal1->GetXaxis()->SetTitle("SigLinPerSigCor");
02176     hSigLinPerSigCorReCal1->GetXaxis()->CenterTitle();
02177     //hSigLinPerSigCorReCal1->SetBit(TH1::kCanRebin);
02178     
02179     hSigCorPerSigMapReCal1=new TH1F("hSigCorPerSigMapReCal1",
02180                                  "hSigCorPerSigMapReCal1",
02181                                  10000,-1,100);
02182     hSigCorPerSigMapReCal1->SetFillColor(0);
02183     hSigCorPerSigMapReCal1->SetTitle("SigCorPerSigMap");
02184     hSigCorPerSigMapReCal1->GetXaxis()->SetTitle("SigCorPerSigMap");
02185     hSigCorPerSigMapReCal1->GetXaxis()->CenterTitle();
02186     //hSigCorPerSigMapReCal1->SetBit(TH1::kCanRebin);
02187     
02188     hSigMapPerMipReCal1=new TH1F("hSigMapPerMipReCal1",
02189                               "hSigMapPerMipReCal1",
02190                               10100,-10,1000);
02191     hSigMapPerMipReCal1->SetFillColor(0);
02192     hSigMapPerMipReCal1->SetTitle("SigMapPerMip");
02193     hSigMapPerMipReCal1->GetXaxis()->SetTitle("SigMapPerMip");
02194     hSigMapPerMipReCal1->GetXaxis()->CenterTitle();
02195     //hSigMapPerMipReCal1->SetBit(TH1::kCanRebin);
02196     
02197     hMipPerGeVReCal1=new TH1F("hMipPerGeVReCal1","hMipPerGeVReCal1",
02198                            10100,-10,1000);
02199     hMipPerGeVReCal1->SetFillColor(0);
02200     hMipPerGeVReCal1->SetTitle("MipPerGeV");
02201     hMipPerGeVReCal1->GetXaxis()->SetTitle("MipPerGeV");
02202     hMipPerGeVReCal1->GetXaxis()->CenterTitle();
02203     //hMipPerGeVReCal1->SetBit(TH1::kCanRebin);
02204 
02205     hAdcPerPeReCal1=new TH1F("hAdcPerPeReCal1","hAdcPerPeReCal1",
02206                            8080,-20,2000);
02207     hAdcPerPeReCal1->SetFillColor(0);
02208     hAdcPerPeReCal1->SetTitle("AdcPerPe");
02209     hAdcPerPeReCal1->GetXaxis()->SetTitle("AdcPerPe");
02210     hAdcPerPeReCal1->GetXaxis()->CenterTitle();
02211     //hAdcPerPeReCal1->SetBit(TH1::kCanRebin);
02212 
02213     hAdcPerSigDrfReCal1=new TH1F("hAdcPerSigDrfReCal1",
02214                                  "hAdcPerSigDrfReCal1",
02215                                  1100,-1,10);
02216     hAdcPerSigDrfReCal1->SetFillColor(0);
02217     hAdcPerSigDrfReCal1->SetTitle("AdcPerSigDrf");
02218     hAdcPerSigDrfReCal1->GetXaxis()->SetTitle("AdcPerSigDrf");
02219     hAdcPerSigDrfReCal1->GetXaxis()->CenterTitle();
02220     //hAdcPerSigDrfReCal1->SetBit(TH1::kCanRebin);
02221 
02222     hAdcPerSigLinOnlyReCal1=new TH1F("hAdcPerSigLinOnlyReCal1",
02223                                      "hAdcPerSigLinOnlyReCal1",
02224                                      2100,-1,20);
02225     hAdcPerSigLinOnlyReCal1->SetFillColor(0);
02226     hAdcPerSigLinOnlyReCal1->SetTitle("AdcPerSigLinOnly");
02227     hAdcPerSigLinOnlyReCal1->GetXaxis()->SetTitle("AdcPerSigLinOnly");
02228     hAdcPerSigLinOnlyReCal1->GetXaxis()->CenterTitle();
02229     //hAdcPerSigLinOnlyReCal1->SetBit(TH1::kCanRebin);
02230 
02231 
02232 
02233     //stripend 2
02234     hAdcPerSigLinReCal2=new TH1F("hAdcPerSigLinReCal2","hAdcPerSigLinReCal2",
02235                               1100,-1,10);
02236     hAdcPerSigLinReCal2->SetFillColor(0);
02237     hAdcPerSigLinReCal2->SetTitle("AdcPerSigLin");
02238     hAdcPerSigLinReCal2->GetXaxis()->SetTitle("AdcPerSigLin");
02239     hAdcPerSigLinReCal2->GetXaxis()->CenterTitle();
02240     //hAdcPerSigLinReCal2->SetBit(TH1::kCanRebin);
02241     
02242     hSigLinPerSigCorReCal2=new TH1F("hSigLinPerSigCorReCal2",
02243                                  "hSigLinPerSigCorReCal2",
02244                                  1100,-1,10);
02245     hSigLinPerSigCorReCal2->SetFillColor(0);
02246     hSigLinPerSigCorReCal2->SetTitle("SigLinPerSigCor");
02247     hSigLinPerSigCorReCal2->GetXaxis()->SetTitle("SigLinPerSigCor");
02248     hSigLinPerSigCorReCal2->GetXaxis()->CenterTitle();
02249     //hSigLinPerSigCorReCal2->SetBit(TH1::kCanRebin);
02250     
02251     hSigCorPerSigMapReCal2=new TH1F("hSigCorPerSigMapReCal2",
02252                                  "hSigCorPerSigMapReCal2",
02253                                  10000,-1,100);
02254     hSigCorPerSigMapReCal2->SetFillColor(0);
02255     hSigCorPerSigMapReCal2->SetTitle("SigCorPerSigMap");
02256     hSigCorPerSigMapReCal2->GetXaxis()->SetTitle("SigCorPerSigMap");
02257     hSigCorPerSigMapReCal2->GetXaxis()->CenterTitle();
02258     //hSigCorPerSigMapReCal2->SetBit(TH1::kCanRebin);
02259     
02260     hSigMapPerMipReCal2=new TH1F("hSigMapPerMipReCal2",
02261                               "hSigMapPerMipReCal2",
02262                               10100,-10,1000);
02263     hSigMapPerMipReCal2->SetFillColor(0);
02264     hSigMapPerMipReCal2->SetTitle("SigMapPerMip");
02265     hSigMapPerMipReCal2->GetXaxis()->SetTitle("SigMapPerMip");
02266     hSigMapPerMipReCal2->GetXaxis()->CenterTitle();
02267     //hSigMapPerMipReCal2->SetBit(TH1::kCanRebin);
02268     
02269     hMipPerGeVReCal2=new TH1F("hMipPerGeVReCal2","hMipPerGeVReCal2",
02270                            10100,-10,1000);
02271     hMipPerGeVReCal2->SetFillColor(0);
02272     hMipPerGeVReCal2->SetTitle("MipPerGeV");
02273     hMipPerGeVReCal2->GetXaxis()->SetTitle("MipPerGeV");
02274     hMipPerGeVReCal2->GetXaxis()->CenterTitle();
02275     //hMipPerGeVReCal2->SetBit(TH1::kCanRebin);
02276 
02277     hAdcPerPeReCal2=new TH1F("hAdcPerPeReCal2","hAdcPerPeReCal2",
02278                            8080,-20,2000);
02279     hAdcPerPeReCal2->SetFillColor(0);
02280     hAdcPerPeReCal2->SetTitle("AdcPerPe");
02281     hAdcPerPeReCal2->GetXaxis()->SetTitle("AdcPerPe");
02282     hAdcPerPeReCal2->GetXaxis()->CenterTitle();
02283     //hAdcPerPeReCal2->SetBit(TH1::kCanRebin);
02284 
02285     hAdcPerSigDrfReCal2=new TH1F("hAdcPerSigDrfReCal2",
02286                                  "hAdcPerSigDrfReCal2",
02287                                  1100,-1,10);
02288     hAdcPerSigDrfReCal2->SetFillColor(0);
02289     hAdcPerSigDrfReCal2->SetTitle("AdcPerSigLinDrf");
02290     hAdcPerSigDrfReCal2->GetXaxis()->SetTitle("AdcPerSigLinDrf");
02291     hAdcPerSigDrfReCal2->GetXaxis()->CenterTitle();
02292     //hAdcPerSigDrfReCal2->SetBit(TH1::kCanRebin);
02293     
02294     hAdcPerSigLinOnlyReCal2=new TH1F("hAdcPerSigLinOnlyReCal2",
02295                                      "hAdcPerSigLinOnlyReCal2",
02296                                      2100,-1,20);
02297     hAdcPerSigLinOnlyReCal2->SetFillColor(0);
02298     hAdcPerSigLinOnlyReCal2->SetTitle("AdcPerSigLinOnly");
02299     hAdcPerSigLinOnlyReCal2->GetXaxis()->SetTitle("AdcPerSigLinOnly");
02300     hAdcPerSigLinOnlyReCal2->GetXaxis()->CenterTitle();
02301     //hAdcPerSigLinOnlyReCal2->SetBit(TH1::kCanRebin);
02302   }
02303 
02304   Calibrator& cal=Calibrator::Instance();
02305   //VldTimeStamp (const time_t &t, const Int_t nsec)
02306   //VldTimeStamp vldts;
02307   VldTimeStamp vldts(s.TimeSec,0);
02308   VldContext vc(static_cast<Detector::Detector_t>(s.Detector),
02309                 static_cast<SimFlag::ESimFlag>(s.SimFlag),vldts);
02310   cal.ReInitialise(vc);//same as reset
02311   //vc.Print();
02312   //cal.PrintConfig(cout);
02313 
02314   MAXMSG("MeuAnalysis",logLevel,1000)
02315     <<"Recalibrate: track "<<s.VtxPlane<<"->"<<s.EndPlane<<endl; 
02316   
02317   //zero all the current values
02318   for (map<Int_t,MeuHitInfo>::iterator it=plInfo.begin();
02319        it!=plInfo.end();++it){
02320     //cache the current cp
02321     MeuHitInfo& cp=it->second;
02322     
02323     MAXMSG("MeuAnalysis",logLevel,1000)
02324       <<"pl="<<cp.Plane      
02325       <<", st="<<cp.Strip
02326       <<", sigCor="<<cp.SigCor
02327       <<", sigMap="<<cp.SigMap
02328       <<", lpos="<<cp.LPos<<endl;
02329 
02330     cp.SigMap=0;
02331     cp.SigMap1=0;
02332     cp.SigMap2=0;
02333     cp.MCSigMap=0;
02334 
02335     if (fRecalibrateSigCor) {
02336       MAXMSG("MeuAnalysis",Msg::kInfo,1)
02337         <<"Reapplying strip-to-strip calibration"
02338         <<", s.SimFlag="<<s.SimFlag<<endl;
02339       cp.SigCor=0;
02340       cp.SigCor1=0;
02341       cp.SigCor2=0;
02342       cp.SigCorTrk1=0;
02343       cp.SigCorTrk2=0;
02344     }
02345     else {
02346       MAXMSG("MeuAnalysis",Msg::kInfo,1)
02347         <<"NOT Reapplying strip-to-strip calibration"<<endl;
02348     }
02349 
02350     if (fRecalibrateSigLin) {
02351       MAXMSG("MeuAnalysis",Msg::kInfo,1)
02352         <<"Reapplying SigLin calibration"<<endl;
02353       cp.SigLin=0;
02354       cp.SigLin1=0;
02355       cp.SigLin2=0;
02356     }
02357     else {
02358       MAXMSG("MeuAnalysis",Msg::kInfo,1)
02359         <<"NOT Reapplying SigLin calibration (hence no drift)"<<endl;
02360     }
02361 
02362     if (fRecalibratePE) {
02363       MAXMSG("MeuAnalysis",Msg::kInfo,1)
02364         <<"Reapplying PE calibration"<<endl;
02365       cp.Pe=0;
02366       cp.Pe1=0;
02367       cp.Pe2=0;
02368     }
02369     else {
02370       MAXMSG("MeuAnalysis",Msg::kInfo,1)
02371         <<"NOT Reapplying PE calibration"<<endl;
02372     }
02373 
02374     //these should already be zero... but to be on the safe side:
02375     cp.SigLinOnly=0;
02376     cp.SigLinOnly1=0;
02377     cp.SigLinOnly2=0;
02378     cp.SigDrf=0;
02379     cp.SigDrf1=0;
02380     cp.SigDrf2=0;
02381   }
02382   
02383   TClonesArray& slcTca=(*ntp.slc);
02384   TClonesArray& stpTca=(*ntp.stp);
02385 
02386   //get the slice associated with this event
02387   const NtpSRSlice* pslc=
02388     dynamic_cast<NtpSRSlice*>(slcTca[evt.slc]);
02389   const NtpSRSlice& slc=(*pslc);
02390   
02392   //loop over strips in slc
02394   for (Int_t i=0;i<slc.nstrip;++i) {
02395     const NtpSRStrip* pstp=
02396       dynamic_cast<NtpSRStrip*>(stpTca[slc.stp[i]]);
02397     const NtpSRStrip& stp=(*pstp);
02398     
02399     Int_t pl=stp.plane;
02400     //check if the plane is in the track
02401     //this cut is not in NuCuts::ExtractPlInfo but it should
02402     //never matter since the track window is always within the track
02403     if (!this->PlaneIsWithinTrack(s,pl,
02404                                   s.VtxPlane,s.EndPlane)) continue;
02405 
02406     //filter out the times that are way off to make the 
02407     //Near and Far detectors more similar
02408     //get the time0, if between 0 and 1 switch to time1
02409     Double_t time=stp.time0;
02410     if (time<0 || time>1) time=stp.time1;
02411     if (time<0 || time>1) cout<<"Ahhhhh, time"<<endl;;
02412     
02413     if (time>s.MedianTime+(400*Munits::ns) ||
02414         time<s.MedianTime-(400*Munits::ns)) {
02415       MAXMSG("MeuCuts",Msg::kInfo,10)
02416         <<"Cut strip: time="<<time
02417         <<", pl="<<stp.plane<<", st="<<stp.strip
02418         <<", sigCor="<<stp.ph0.sigcor+stp.ph1.sigcor
02419         <<", dT="<<(time-s.MedianTime)/Munits::ns<<endl;
02420       continue;//messes up beg/end planes (if not filtered bad trks)!
02421     }
02422 
02423     //cache the current cp
02424     MeuHitInfo& cp=plInfo[pl];
02425 
02426     //get the raw adcs
02427     Float_t adc1=stp.ph0.raw;//east
02428     Float_t adc2=stp.ph1.raw;//west
02429 
02430     PlexStripEndId seidE(static_cast<Detector::Detector_t>
02431                          (s.Detector),pl,stp.strip,
02432                          StripEnd::kEast);
02433     PlexStripEndId seidW(static_cast<Detector::Detector_t>
02434                          (s.Detector),pl,stp.strip,
02435                          StripEnd::kWest);
02436     
02437     //sigDrf and sigLinOnly values are NOT used to get the MEU number
02438     //get the sigLinOnly value (applied to strip)
02439     //useful for systematic error evaluation
02440     //NOTE: in macro have to set:
02441     //cal.LinearityCalScheme().Set("BucketMode=0"); 
02442     //for the results to have any meaning
02443     Float_t sigLinOnly1=0;
02444     Float_t sigLinOnly2=0;
02445     if (adc1>0) sigLinOnly1=cal.GetLinearized(adc1,seidE);
02446     else sigLinOnly1=0;
02447     if (adc2>0) sigLinOnly2=cal.GetLinearized(adc2,seidW);
02448     else sigLinOnly2=0;
02449 
02450     //now sum up everything in the plane
02451     cp.SigLinOnly+=sigLinOnly1+sigLinOnly2;
02452     cp.SigLinOnly1+=sigLinOnly1;
02453     cp.SigLinOnly2+=sigLinOnly2;
02454     
02455     //get the drift correction applied to adcs
02456     Float_t sigDrf1=0;
02457     Float_t sigDrf2=0;
02458     if (adc1>0) sigDrf1=cal.GetDriftCorrected(adc1,seidE);
02459     else sigDrf1=0;
02460     if (adc2>0) sigDrf2=cal.GetDriftCorrected(adc2,seidW);
02461     else sigDrf2=0;
02462 
02463     //now sum up everything in the plane
02464     cp.SigDrf+=sigDrf1+sigDrf2;
02465     cp.SigDrf1+=sigDrf1;
02466     cp.SigDrf2+=sigDrf2;
02467 
02468 
02470     //section to recalibrate PEs if required
02471     Float_t pe1=stp.ph0.pe;//east
02472     Float_t pe2=stp.ph1.pe;//west
02473     if (fRecalibratePE){
02474       if (adc1>0) pe1=cal.GetPhotoElectrons(adc1,seidE);
02475       else pe1=0;
02476       if (adc2>0) pe2=cal.GetPhotoElectrons(adc2,seidW);
02477       else pe2=0;
02478 
02479       //now sum up everything in the plane
02480       cp.Pe+=pe1+pe2;
02481       cp.Pe1+=pe1;
02482       cp.Pe2+=pe2;
02483     }
02484     
02486     //section to recalibrate siglin if required
02487     Float_t sigLin1=stp.ph0.siglin;//east
02488     Float_t sigLin2=stp.ph1.siglin;//west
02489     if (fRecalibrateSigLin){
02490       //linearity calibration is done on digits (not strips)
02491       //this matters for the ND!
02492       //It is not totally correct to recalibrate with the standard
02493       //ntuples because the digit information is not available for ND
02494       //however the calibrator can implement an accurate fudge:
02495       //cal.LinearityCalScheme().Set("BucketMode=0"); 
02496 
02497       MAXMSG("MeuAnalysis",Msg::kInfo,1)
02498         <<"Recalibrating drift and linearity for detector="
02499         <<s.Detector<<endl;
02500       
02501       //first apply the linearity 
02502       if (adc1>0) sigLin1=cal.GetLinearized(adc1,seidE);
02503       else sigLin1=0;
02504       if (adc2>0) sigLin2=cal.GetLinearized(adc2,seidW);
02505       else sigLin2=0;
02506       
02507       //then do the drift correction (to sigLin - driftCor this time)
02508       if (sigLin1>0) sigLin1=cal.GetDriftCorrected(sigLin1,seidE);
02509       else sigLin1=0;
02510       if (sigLin2>0) sigLin2=cal.GetDriftCorrected(sigLin2,seidW);
02511       else sigLin2=0;   
02512       
02513       //now sum up everything in the plane
02514       cp.SigLin+=sigLin1+sigLin2;
02515       cp.SigLin1+=sigLin1;
02516       cp.SigLin2+=sigLin2;
02517     }
02518     
02520     //section to recalibrate sigcor if required
02521     Float_t sigCor1=stp.ph0.sigcor;//east
02522     Float_t sigCor2=stp.ph1.sigcor;//west
02523     if (fRecalibrateSigCor){
02524       MAXMSG("MeuAnalysis",Msg::kInfo,1)
02525         <<"Recalibrating strip-to-strip..."<<endl;
02526       if (s.Detector!=Detector::kNear) {
02527         if (sigLin1>0) sigCor1=cal.GetStripToStripCorrected(sigLin1,
02528                                                             seidE);
02529       }
02530       //else sigCor1=0;
02531       if (sigLin2>0) sigCor2=cal.GetStripToStripCorrected(sigLin2,
02532                                                           seidW);
02533       else sigCor2=0;
02534       
02535       //sum both stripends
02536       Float_t sigCorBoth=sigCor1+sigCor2;
02537       
02538       //sum the sigcor in the object
02539       cp.SigCor+=sigCorBoth;
02540       cp.SigCor1+=sigCor1;
02541       cp.SigCor2+=sigCor2;
02542       
02543       //store the sigcor of both ends of the tracked strip
02544       //this means that fSigCor!=fSigCor1+fSigCor2 always
02545       if (stp.strip==cp.Strip){
02546         cp.SigCorTrk1=sigCor1;
02547         cp.SigCorTrk2=sigCor2;
02548       }
02549     }
02550     else if (fRecalibrateSigLin) {
02551       MAXMSG("MeuAnalysis",Msg::kInfo,1)
02552         <<"Using ntp strip-to-strip calibration on recalibrated"
02553         <<" SigLin values..."<<endl;
02554       //if siglin was recalibrated then have to feed the value through
02555       //the value of the sigcor correction in the ntp is just
02556       //stp.ph0.sigcor/stp.ph0.siglin so multiply by that:
02557       sigCor1=sigLin1*stp.ph0.sigcor/stp.ph0.siglin;//east
02558       sigCor2=sigLin2*stp.ph1.sigcor/stp.ph1.siglin;//west
02559     }
02560 
02561     //find the lpos for this plane then calibrate
02562     Float_t lpos=cp.LPos;
02563     Float_t sigMap1=0;
02564     Float_t sigMap2=0;
02565 
02566     //only calibrate the east side for !ND
02567     //also only do the calibration if there is a hit on stripend
02568     if (s.Detector!=Detector::kNear) {
02569       if (sigCor1>0) sigMap1=cal.GetAttenCorrectedTpos(sigCor1,
02570                                                        lpos,seidE);
02571     }
02572     if (sigCor2>0) sigMap2=cal.GetAttenCorrectedTpos(sigCor2,
02573                                                      lpos,seidW);
02574     
02575     Float_t sigMapBoth=sigMap1+sigMap2;
02576     
02577     //sum the sigmap
02578     cp.SigMap+=sigMapBoth;
02579     cp.SigMap1+=sigMap1;
02580     cp.SigMap2+=sigMap2;
02581     
02582     MAXMSG("MeuAnalysis",logLevel,1000)
02583       <<"strip="<<stp.strip<<", pl="<<pl<<", lpos="<<lpos<<endl
02584       <<"  sigCor1="<<sigCor1<<", sigMap1="<<sigMap1<<endl
02585       <<"  sigCor2="<<sigCor2<<", sigMap2="<<sigMap2<<endl
02586       <<endl;
02587 
02588     //calibrate "true MC" values as well
02589     if (s.SimFlag==SimFlag::kMC){
02590       //get the MC lpos
02591       Float_t lpos=cp.MCLPos;
02592 
02593       Float_t sigMap1=0;
02594       Float_t sigMap2=0;
02595 
02596       if (sigCor1>0) {
02597         sigMap1=cal.GetAttenCorrectedTpos(sigCor1,lpos,seidE);
02598       }
02599       if (sigCor2>0) {
02600         sigMap2=cal.GetAttenCorrectedTpos(sigCor2,lpos,seidW);
02601       }
02602 
02603       //calc sigmap and add to cp
02604       sigMapBoth=sigMap1+sigMap2;
02605       cp.MCSigMap+=sigMapBoth;
02606     }
02607 
02608     //fill calibration value histos
02609     //can't do mip or gev ones
02610     if (sigLin1) hAdcPerSigLinReCal1->Fill(adc1/sigLin1);
02611     if (sigCor1) hSigLinPerSigCorReCal1->Fill(sigLin1/sigCor1);
02612     if (sigCor1) hSigCorPerSigMapReCal1->Fill(sigCor1/sigMap1);
02613     if (pe1) hAdcPerPeReCal1->Fill(adc1/pe1);
02614     if (sigLinOnly1) hAdcPerSigLinOnlyReCal1->Fill(adc1/sigLinOnly1);
02615     if (sigDrf1) hAdcPerSigDrfReCal1->Fill(adc1/sigDrf1);
02616     //stripend 2
02617     if (sigLin2) hAdcPerSigLinReCal2->Fill(adc2/sigLin2);
02618     if (sigCor2) hSigLinPerSigCorReCal2->Fill(sigLin2/sigCor2);
02619     if (sigCor2) hSigCorPerSigMapReCal2->Fill(sigCor2/sigMap2);
02620     if (pe2) hAdcPerPeReCal2->Fill(adc2/pe2);
02621     if (sigLinOnly2) hAdcPerSigLinOnlyReCal2->Fill(adc2/sigLinOnly2);
02622     if (sigDrf2) hAdcPerSigDrfReCal2->Fill(adc2/sigDrf2);
02623     
02624     //print out
02625     static Bool_t print=false;
02626     if (print) {
02627       if (sigDrf2) {
02628         if (adc2/sigDrf2>2) {
02629           MAXMSG("MeuAnalysis",Msg::kInfo,1000)
02630             <<"adc2/sigDrf2="<<adc2/sigDrf2
02631             <<", adc2="<<adc2<<", sigDrf2="<<sigDrf2<<endl;
02632         }
02633       }
02634       if (sigLinOnly2) {
02635         if (adc2/sigLinOnly2>2) {
02636           MAXMSG("MeuAnalysis",Msg::kInfo,1000)
02637             <<"adc2/sigLinOnly2="<<adc2/sigLinOnly2
02638             <<", adc2="<<adc2<<", sigLinOnly2="<<sigLinOnly2
02639             <<", pl,st="<<cp.Plane<<","<<cp.Strip<<endl;
02640         }
02641       }
02642     }
02643   }//loop over strips in slc
02644 
02645   //print
02646   Bool_t printMap=false;
02647   if (printMap){
02648     for (map<Int_t,MeuHitInfo>::iterator it=plInfo.begin();
02649          it!=plInfo.end();++it){
02650       //cache the current cp
02651       MeuHitInfo& cp=it->second;
02652       
02653       MAXMSG("MeuAnalysis",Msg::kInfo,1000)
02654         <<"Printing: pl="<<cp.Plane      
02655         <<", st="<<cp.Strip
02656         <<", sigCor="<<cp.SigCor
02657         <<", sigMap="<<cp.SigMap
02658         <<", lpos="<<cp.LPos<<endl;
02659     }
02660   }
02661 }

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

Definition at line 2666 of file MeuAnalysis.cxx.

References MeuSummary::Detector, MeuSummary::EndPlane, Calibrator::GetAttenCorrectedTpos(), Calibrator::GetDriftCorrected(), Calibrator::GetLinearized(), Calibrator::GetPhotoElectrons(), Calibrator::GetStripToStripCorrected(), Calibrator::Instance(), it, Msg::kDebug, StripEnd::kEast, Detector::kFar, Msg::kInfo, SimFlag::kMC, Detector::kNear, StripEnd::kWest, 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(), CfgPromptConfigurable::Set(), MeuHitInfo::SigCor, MeuHitInfo::SigCor1, MeuHitInfo::SigCor2, MeuHitInfo::SigCorTrk1, MeuHitInfo::SigCorTrk2, MeuHitInfo::SigLin, MeuHitInfo::SigLin1, MeuHitInfo::SigLin2, MeuHitInfo::SigMap, MeuHitInfo::SigMap1, MeuHitInfo::SigMap2, MeuSummary::SimFlag, AtmosStrip::Strip, MeuHitInfo::Strip, AtmosEvent::StripList, AtmosEvent::UnixTime, and MeuSummary::VtxPlane.

Referenced by MakeSummaryTreeWithAtNu().

02668 {
02669   //variable to turn on all the useful messages if required
02670   Msg::LogLevel_t logLevel=Msg::kDebug;
02671   //Msg::LogLevel_t logLevelSts=Msg::kDebug;
02672   //Msg::LogLevel_t logLevelLin=Msg::kInfo;
02673 
02674   Bool_t fRecalibrateSigLin=true;//have to configure the calibrator too! 
02675   Bool_t fRecalibrateSigCor=true;//have to configure the calibrator too! 
02676   Bool_t fRecalibratePE=true;//have to configure the calibrator too! 
02677 
02678   Calibrator& cal=Calibrator::Instance();
02679   //VldTimeStamp (const time_t &t, const Int_t nsec)
02680   //VldTimeStamp vldts;
02681   VldTimeStamp vldts(ev.UnixTime,0);
02682   VldContext vc(static_cast<Detector::Detector_t>(s.Detector),
02683                 static_cast<SimFlag::ESimFlag>(s.SimFlag),vldts);
02684   cal.ReInitialise(vc);
02685 
02686   static Bool_t firstTime=false;//switch this off with false
02687   if (firstTime){
02688     if (s.Detector==Detector::kFar){
02689       cal.Set("LinCalibrator=PulserLinearityCalScheme");
02690     }
02691     else if (s.Detector==Detector::kNear){
02692       cal.Set("LinCalibrator=QuadLinearityCalScheme"); 
02693     }
02694     else cout<<"Ahh, detector="<<s.Detector<<endl;
02695     
02696     cout<<"Actually using calibrator configured as:"<<endl;
02697     cal.PrintConfig(cout);
02698     firstTime=false;
02699   }
02700   
02701   MAXMSG("MeuAnalysis",logLevel,1000)
02702     <<"Recalibrate: track "<<s.VtxPlane<<"->"<<s.EndPlane<<endl; 
02703 
02704   //zero all the current values
02705   for (map<Int_t,MeuHitInfo>::iterator it=plInfo.begin();
02706        it!=plInfo.end();++it){
02707     //cache the current cp
02708     MeuHitInfo& cp=it->second;
02709     
02710     MAXMSG("MeuAnalysis",logLevel,1000)
02711       <<"pl="<<cp.Plane      
02712       <<", st="<<cp.Strip
02713       <<", sigCor="<<cp.SigCor
02714       <<", sigMap="<<cp.SigMap
02715       <<", lpos="<<cp.LPos<<endl;
02716 
02717     cp.SigMap=0;
02718     cp.SigMap1=0;
02719     cp.SigMap2=0;
02720     cp.MCSigMap=0;
02721     if (fRecalibrateSigCor) {
02722       MAXMSG("MeuAnalysis",Msg::kInfo,1)
02723         <<"Reapplying strip-to-strip calibration, cp.SimFlag="
02724         <<s.SimFlag<<endl;
02725       cp.SigCor=0;
02726       cp.SigCor1=0;
02727       cp.SigCor2=0;
02728       cp.SigCorTrk1=0;
02729       cp.SigCorTrk2=0;
02730     }
02731     else {
02732       MAXMSG("MeuAnalysis",Msg::kInfo,1)
02733         <<"NOT Reapplying strip-to-strip calibration"<<endl;
02734     }
02735 
02736     if (fRecalibrateSigLin) {
02737       MAXMSG("MeuAnalysis",Msg::kInfo,1)
02738         <<"Reapplying SigLin calibration"<<endl;
02739       cp.SigLin=0;
02740       cp.SigLin1=0;
02741       cp.SigLin2=0;
02742     }
02743     else {
02744       MAXMSG("MeuAnalysis",Msg::kInfo,1)
02745         <<"NOT Reapplying SigLin calibration"<<endl;
02746     }
02747 
02748     if (fRecalibratePE) {
02749       MAXMSG("MeuAnalysis",Msg::kInfo,1)
02750         <<"Reapplying PE calibration"<<endl;
02751       cp.Pe=0;
02752       cp.Pe1=0;
02753       cp.Pe2=0;
02754     }
02755     else {
02756       MAXMSG("MeuAnalysis",Msg::kInfo,1)
02757         <<"NOT Reapplying PE calibration"<<endl;
02758     }
02759 
02760   }
02761 
02762   TClonesArray& strips=(*ev.StripList);
02763   Int_t numStrips=strips.GetEntries();
02765   //loop over the strips in the snarl
02767   for (Int_t hit=0;hit<numStrips;hit++){
02768     const AtmosStrip* strip=dynamic_cast<AtmosStrip*>(strips[hit]);
02769     const AtmosStrip& stp=(*strip);
02770     
02771     Int_t pl=stp.Plane;
02772     //check if the plane is in the track
02773     if (!this->PlaneIsWithinTrack(s,pl,
02774                                   s.VtxPlane,s.EndPlane)) continue;
02775 
02776     //cache the current cp
02777     MeuHitInfo& cp=plInfo[pl];
02778 
02779     Float_t adc1=stp.Qadc[0];//east
02780     Float_t adc2=stp.Qadc[1];//west
02781     Float_t pe1=stp.QPE[0];//east
02782     Float_t pe2=stp.QPE[1];//west
02783     Float_t sigLin1=stp.Qadc[0];//east//siglin is no longer in ntp
02784     Float_t sigLin2=stp.Qadc[1];//west//so use adcs instead
02785     Float_t sigCor1=stp.Qadc[0];//east//siglin is no longer in ntp
02786     Float_t sigCor2=stp.Qadc[1];//west//so use adcs instead
02787     //Float_t sigCor1=stp.QPEcorr[0];//east//this is pe size so don't use it!
02788     //Float_t sigCor2=stp.QPEcorr[1];//west
02789     
02790     PlexStripEndId seidE(static_cast<Detector::Detector_t>
02791                          (s.Detector),pl,stp.Strip,
02792                          StripEnd::kEast);
02793     PlexStripEndId seidW(static_cast<Detector::Detector_t>
02794                          (s.Detector),pl,stp.Strip,
02795                          StripEnd::kWest);
02796 
02798     //section to recalibrate PEs if required
02799     if (fRecalibratePE){
02800       if (adc1>0) pe1=cal.GetPhotoElectrons(adc1,seidE);
02801       else pe1=0;
02802       if (adc2>0) pe2=cal.GetPhotoElectrons(adc2,seidW);
02803       else pe2=0;
02804 
02805       //now sum up everything in the plane
02806       cp.Pe+=pe1+pe2;
02807       cp.Pe1+=pe1;
02808       cp.Pe2+=pe2;
02809     }
02810 
02812     //section to recalibrate siglin if required
02813     if (fRecalibrateSigLin){
02814       //first do the drift correction
02815       if (adc1>0) sigLin1=cal.GetDriftCorrected(adc1,seidE);
02816       else sigLin1=0;
02817       if (adc2>0) sigLin2=cal.GetDriftCorrected(adc2,seidW);
02818       else sigLin2=0;
02819 
02820       //then apply the linearity (to sigLin - driftCor this time)
02821       if (adc1>0) sigLin1=cal.GetLinearized(sigLin1,seidE);
02822       else sigLin1=0;
02823       if (adc2>0) sigLin2=cal.GetLinearized(sigLin2,seidW);
02824       else sigLin2=0;
02825 
02826       //now sum up everything in the plane
02827       cp.SigLin+=sigLin1+sigLin2;
02828       cp.SigLin1+=sigLin1;
02829       cp.SigLin2+=sigLin2;
02830     }
02831 
02833     //section to recalibrate sigcor if required
02834     //  have to do this for AtNu for MC AND DATA 
02835     //  since it has pe-sized units of "cor"
02836     if (fRecalibrateSigCor){// && s.SimFlag==SimFlag::kData){
02837       if (sigLin1>0) sigCor1=cal.GetStripToStripCorrected(sigLin1,
02838                                                           seidE);
02839       else sigCor1=0;
02840       if (sigLin2>0) sigCor2=cal.GetStripToStripCorrected(sigLin2,
02841                                                           seidW);
02842       else sigCor2=0;
02843       
02844       //sum the sigcor in the object
02845       cp.SigCor+=sigCor1+sigCor2;
02846       cp.SigCor1+=sigCor1;
02847       cp.SigCor2+=sigCor2;
02848       
02849       //store the sigcor of both ends of the tracked strip
02850       //this means that fSigCor!=fSigCor1+fSigCor2 always
02851       if (stp.Strip==cp.Strip){
02852         cp.SigCorTrk1=sigCor1;
02853         cp.SigCorTrk2=sigCor2;
02854       }
02855     }
02856 
02857     //find the lpos for this plane then calibrate
02858     Float_t lpos=cp.LPos;
02859     Float_t sigMap1=0;
02860     Float_t sigMap2=0;
02861 
02862     //only do the look up if non zero
02863     if (sigCor1>0) sigMap1=cal.GetAttenCorrectedTpos(sigCor1,
02864                                                      lpos,seidE);
02865     if (sigCor2>0) sigMap2=cal.GetAttenCorrectedTpos(sigCor2,
02866                                                      lpos,seidW);
02867     
02868     Float_t sigMapBoth=sigMap1+sigMap2;
02869     
02870     //sum the sigmap
02871     cp.SigMap+=sigMapBoth;
02872     //store the sigmap of both ends of the tracked strip
02873     //this means that fSigMap!=fSigMap1+fSigMap2 always
02874     //if (st.Strip==cp.Strip){
02875     cp.SigMap1+=sigMap1;
02876     cp.SigMap2+=sigMap2;
02877     //}
02878 
02879     MAXMSG("MeuAnalysis",logLevel,1000)
02880       <<"strip="<<stp.Strip<<", pl="<<pl<<", lpos="<<lpos<<endl
02881       <<"  sigCor1="<<sigCor1<<", sigMap1="<<sigMap1<<endl
02882       <<"  sigCor2="<<sigCor2<<", sigMap2="<<sigMap2<<endl
02883       <<endl;
02884 
02885     //do MC as well
02886     if (s.SimFlag==SimFlag::kMC){
02887       //get the MC lpos
02888       Float_t lpos=cp.MCLPos;
02889 
02890       Float_t sigMap1=0;
02891       Float_t sigMap2=0;
02892 
02893       if (sigCor1>0) {
02894         sigMap1=cal.GetAttenCorrectedTpos(sigCor1,lpos,seidE);
02895       }
02896       if (sigCor2>0) {
02897         sigMap2=cal.GetAttenCorrectedTpos(sigCor2,lpos,seidW);
02898       }
02899 
02900       //calc sigmap and add to cp
02901       sigMapBoth=sigMap1+sigMap2;
02902       cp.MCSigMap+=sigMapBoth;
02903     }
02904   }
02905 
02906   //print
02907   Bool_t printMap=false;
02908   if (printMap){
02909     for (map<Int_t,MeuHitInfo>::iterator it=plInfo.begin();
02910          it!=plInfo.end();++it){
02911       //cache the current cp
02912       MeuHitInfo& cp=it->second;
02913       
02914       MAXMSG("MeuAnalysis",Msg::kInfo,1000)
02915         <<"Printing: pl="<<cp.Plane      
02916         <<", st="<<cp.Strip
02917         <<", sigCor="<<cp.SigCor
02918         <<", sigMap="<<cp.SigMap
02919         <<", lpos="<<cp.LPos<<endl;
02920     }
02921   }
02922 }

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

Definition at line 201 of file MeuAnalysis.cxx.

References fRecalibratePE, Msg::kInfo, and MSG.

00202 {
00203   MSG("MeuAnalysis",Msg::kInfo)
00204     <<"Setting fRecalibratePE="<<recalibrate<<endl;
00205   fRecalibratePE=recalibrate;
00206   if (recalibrate) {
00207     MSG("MeuAnalysis",Msg::kInfo)
00208       <<"Note that the calibrator has to be correctly configured"
00209       <<" for the recalibration to work!"<<endl;
00210   }
00211 }

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

Definition at line 187 of file MeuAnalysis.cxx.

References fRecalibrateSigCor, Msg::kInfo, and MSG.

00188 {
00189   MSG("MeuAnalysis",Msg::kInfo)
00190     <<"Setting fRecalibrateSigCor="<<recalibrate<<endl;
00191   fRecalibrateSigCor=recalibrate;
00192   if (recalibrate) {
00193     MSG("MeuAnalysis",Msg::kInfo)
00194       <<"Note that the calibrator has to be correctly configured"
00195       <<" for the recalibration to work!"<<endl;
00196   }
00197 }

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

Definition at line 173 of file MeuAnalysis.cxx.

References fRecalibrateSigLin, Msg::kInfo, and MSG.

00174 {
00175   MSG("MeuAnalysis",Msg::kInfo)
00176     <<"Setting fRecalibrateSigLin="<<recalibrate<<endl;
00177   fRecalibrateSigLin=recalibrate;
00178   if (recalibrate) {
00179     MSG("MeuAnalysis",Msg::kInfo)
00180       <<"Note that the calibrator has to be correctly configured"
00181       <<" for the recalibration to work!"<<endl;
00182   }
00183 }

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

Definition at line 743 of file MeuAnalysis.cxx.

References fChain, fChainBD, fEntries, Msg::kInfo, and MSG.

Referenced by BasicPlots(), BasicReco(), MakeSummaryTreeWithAtNu(), MakeSummaryTreeWithNtpSt(), N_1Plots(), SnarlList(), SpillPlots(), and TestEventLoop().

00744 {
00745   if (printMode==1){
00746     Float_t fract=ceil(fEntries/20.);  
00747     if (ceil(((Float_t)entry)/fract)==((Float_t)entry)/fract){
00748       MSG("MeuAnalysis",Msg::kInfo) 
00749         <<"Fraction of loop complete: "<<entry 
00750         <<"/"<<fEntries<<"  ("
00751         <<(Int_t)(100.*entry/fEntries)<<"%)"<<endl;
00752     }
00753   }
00754 
00755   //get the snarl/entry
00756   fChain->GetEvent(entry);
00757   if (fChainBD) fChainBD->GetEvent(entry); 
00758 }

void MeuAnalysis::SetRootFileObjects (  )  [private]

Definition at line 593 of file MeuAnalysis.cxx.

References NtpStRecord::evthdr, fChain, fChainBD, fEntries, fEntriesBD, fev, fRec, fRecBD, fUseAtNu, Msg::kDebug, Msg::kInfo, MSG, NtpSREventSummary::nslice, and NtpBDLiteRecord::tortgt.

00594 {
00595   MSG("MeuAnalysis",Msg::kInfo)
00596     <<"Running the SetRootFileObjects method..."<<endl;
00597 
00598   if (fUseAtNu) {
00599     fChain->SetBranchAddress("evt", &fev);
00600     
00601     //set the number of events in the tree
00602     fEntries=static_cast<Int_t>(fChain->GetEntries());
00603     MSG("MeuAnalysis",Msg::kInfo)
00604       <<"AtNuTree has "<<fEntries<<" entries"<<endl;
00605     if (fEntries>0){
00606       //get the snarl/event
00607       fChain->GetEvent(0); 
00608       //fRunNumber=(*fev).Run;
00609       MSG("MeuAnalysis",Msg::kInfo)
00610         <<"SetRootFileObjects: first UnixTime="<<(*fev).UnixTime<<endl;
00611     }
00612   }
00613   else {
00614     //set up tree
00615     fChain->SetBranchAddress("NtpStRecord",&fRec);
00616     if (fChainBD) fChainBD->SetBranchAddress("NtpBDLiteRecord",&fRecBD);
00617     //fBDtortgt
00618     
00619     //set the number of events in the tree
00620     fEntries=static_cast<Int_t>(fChain->GetEntries());
00621     if (fChainBD) fEntriesBD=static_cast<Int_t>(fChainBD->GetEntries());
00622     MSG("MeuAnalysis",Msg::kInfo)
00623       <<"NtpSt tree has "<<fEntries<<" entries"<<endl;
00624     MSG("MeuAnalysis",Msg::kInfo)
00625       <<"NtpBDLite tree has "<<fEntriesBD<<" entries"<<endl;
00626     
00627     if (fEntries>0){
00628       //get the snarl/event
00629       fChain->GetEvent(0); 
00630       //MSG("MeuAnalysis",Msg::kInfo)
00631       //<<"First run number="<<fRec->fHeader.fRun<<endl;
00632       MSG("MeuAnalysis",Msg::kInfo)
00633         <<"SetRootFileObjects: first snarl nslice="
00634         <<fRec->evthdr.nslice<<endl;
00635     }
00636 
00637     if (fEntriesBD>0){
00638       //get the snarl/event
00639       fChainBD->GetEvent(0); 
00640       MSG("MeuAnalysis",Msg::kInfo)
00641         <<"First tortgt="<<fRecBD->tortgt<<" E12 POT"<<endl;
00642     }
00643   }
00644   MSG("MeuAnalysis",Msg::kDebug)
00645     <<"Finished the SetRootFileObjects method"<<endl;
00646 }

void MeuAnalysis::SnarlList (  ) 

Definition at line 3183 of file MeuAnalysis.cxx.

References MeuReco::CalcFCPC(), MeuReco::CalcVtx(), MeuReco::CalcVtxSpecialND(), MeuSummary::DistToEdgeFid, NtpSRTrack::end, MeuSummary::EndPlane, MeuSummary::Event, NtpStRecord::evt, MeuCuts::ExtractPlInfo(), fEntries, MeuCuts::FillSTSumDetails(), RecRecordImp< T >::GetHeader(), RecDataHeader::GetRun(), RecPhysicsHeader::GetSnarl(), InitialiseLoopVariables(), Msg::kDebug, Msg::kInfo, MAXMSG, MSG, OpenTxtFile(), MeuSummary::PC, NtpSRVertex::plane, MeuSummary::Reset(), MeuSummary::RFid, MeuSummary::Run, SetLoopVariables(), MeuSummary::Snarl, NtpStRecord::trk, and NtpSRTrack::vtx.

03184 {
03185   MSG("MeuAnalysis",Msg::kInfo) 
03186     <<" ** Running SnarlList method... **"<<endl;
03187 
03188   ofstream& snarlList=*(this->OpenTxtFile(100,"snarlList"));
03189   
03190   //reconstruction object
03191   MeuReco reco;
03192   MeuCuts cuts;
03193   MeuSummary s;
03194   
03195   //variable to turn on all the useful messages if required
03196   Msg::LogLevel_t logLevel=Msg::kDebug;
03197 
03201   
03202   this->InitialiseLoopVariables();  
03203   
03204   //for(Int_t entry=55000;entry<fEntries;entry++){
03205   for(Int_t entry=0;entry<fEntries;entry++){
03206       
03207     this->SetLoopVariables(entry);
03208     
03209     MSG("MeuAnalysis",logLevel)
03210       <<"Entry="<<entry<<endl;
03211 
03212     //make a reference instead of a pointer
03213     NtpStRecord& ntp=(*fRec);
03214     const RecCandHeader& rec=ntp.GetHeader();
03215     MAXMSG("MeuAnalysis",Msg::kInfo,1)
03216       <<"First entry: run="<<rec.GetRun()
03217       <<", snarl="<<rec.GetSnarl()<<endl;
03218 
03219     //MeuSummary s;
03220     s.Reset();
03221     s.RFid=1.0;
03222     s.DistToEdgeFid=0.4;
03223     s.Event=entry;
03224     cuts.FillSTSumDetails(ntp,s,-1,-1);
03225     //a map to hold the info for each plane in the entry
03226     map<Int_t,MeuHitInfo> plInfo;
03227 
03228     //ensure there is only 1 track in the entry
03229     TClonesArray& trkTca=(*ntp.trk);
03230     Int_t numTrks=trkTca.GetEntries();
03231     if (numTrks!=1) continue;
03232     const NtpSRTrack* ptrk=dynamic_cast<NtpSRTrack*>(trkTca[0]);
03233     const NtpSRTrack& trk=(*ptrk);
03234     
03235     //cut on number of planes
03236     Int_t numPlanes=abs(trk.vtx.plane-trk.end.plane);
03237     if (numPlanes<6) continue;
03238 
03239     TClonesArray& evtTca=(*ntp.evt);
03240     const NtpSREvent* pevt=dynamic_cast<NtpSREvent*>(evtTca[0]);
03241     cuts.ExtractPlInfo(ntp,s,plInfo,*pevt);
03242     
03243     reco.CalcVtx(s,plInfo);
03244 
03245     //check if track stops in spectrometer
03246     if (s.EndPlane>120) continue;
03247 
03248     reco.CalcFCPC(s,plInfo);
03249 
03250     //recalculate the vtx of the tracks coming in from the spectrometer
03251     reco.CalcVtxSpecialND(s,plInfo);
03252 
03253     if (s.PC) {
03254       MAXMSG("MeuAnalysis",Msg::kInfo,100)
03255         <<"Entry="<<entry<<", run="<<rec.GetRun()
03256         <<", snarl="<<rec.GetSnarl()<<endl;
03257       snarlList<<s.Run<<"\t"<<s.Snarl<<endl;
03258     }
03259   }//end of for                                       
03260   
03264 
03265   MSG("MeuAnalysis",Msg::kInfo)<<"Finished main loop"<<endl;
03266 
03267   MSG("MeuAnalysis",Msg::kInfo) 
03268     <<" ** Finished SnarlList method **"<<endl;
03269 }

void MeuAnalysis::SpillPlots (  ) 

<<PCCounter-bothSMCounter-coilHitCounter-

Definition at line 3273 of file MeuAnalysis.cxx.

References 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, fEntries, MeuCuts::FillEnVsDist(), MeuCuts::FillEnVsDist2(), MeuCuts::FillSTSumDetails(), MeuCuts::FillSTSumRecoDetails(), MeuCuts::FillTimeHistos(), MeuCuts::FilterBadDistEndStrip(), MeuCuts::FilterBadEndGaps(), MeuCuts::FilterBadEvtPerSlc(), MeuCuts::FilterBadTrackEnd(), MeuCuts::FilterBadTrkTimes(), fOutFile, fRec, fRun, RecRecordImp< T >::GetHeader(), RecDataHeader::GetRun(), RecPhysicsHeader::GetSnarl(), Munits::GeV, InitialiseLoopVariables(), it, Msg::kDebug, Msg::kError, Msg::kInfo, SimFlag::kMC, Msg::kVerbose, MeuHitInfo::LPos, MAXMSG, MeuSummary::MCEndPlane, MeuSummary::MCHighEn, MeuSummary::MCLowEn, MeuSummary::MCVtxPlane, MSG, NtpSREvent::nshower, NtpSRTrack::nstrip, NtpSREvent::ntrack, OpenFile(), OpenTxtFile(), MeuSummary::PC, MeuHitInfo::Plane, NtpSRShower::plane, NtpSRTrack::plane, print(), MeuCuts::PrintNtpSt(), Recalibrate(), MeuReco::Reconstruct(), MeuSummary::RFid, MeuSummary::Run, SetLoopVariables(), NtpSREvent::shw, NtpStRecord::shw, MeuHitInfo::SigMap, MeuSummary::SimFlag, NtpSREvent::slc, MeuSummary::Snarl, MeuHitInfo::TPos, NtpStRecord::trk, NtpSREvent::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.

03274 {
03275   MSG("MeuAnalysis",Msg::kInfo) 
03276     <<" ** Running SpillPlots method... **"<<endl;
03277 
03278   //open the output file for the histograms
03279   fOutFile=this->OpenFile(100,"SpillPlots");
03280 
03281   TH1F* hWinSigMap=new TH1F("hWinSigMap","hWinSigMap",3000,-1,4000);
03282   hWinSigMap->SetTitle("hWinSigMap");
03283   hWinSigMap->GetXaxis()->SetTitle("hWinSigMap");
03284   hWinSigMap->GetXaxis()->CenterTitle();
03285   hWinSigMap->SetFillColor(0);
03286   hWinSigMap->SetLineWidth(3);
03287   //hWinSigMap->SetBit(TH1::kCanRebin);
03288 
03289   TH1F* hWinSigCor=new TH1F("hWinSigCor","hWinSigCor",3000,-1,4000);
03290   hWinSigCor->SetTitle("hWinSigCor");
03291   hWinSigCor->GetXaxis()->SetTitle("hWinSigCor");
03292   hWinSigCor->GetXaxis()->CenterTitle();
03293   hWinSigCor->SetFillColor(0);
03294   hWinSigCor->SetLineWidth(3);
03295   //hWinSigCor->SetBit(TH1::kCanRebin);
03296   
03297   TH1F* hWinSigLin=new TH1F("hWinSigLin","hWinSigLin",3000,-1,4000);
03298   hWinSigLin->SetTitle("hWinSigLin");
03299   hWinSigLin->GetXaxis()->SetTitle("hWinSigLin");
03300   hWinSigLin->GetXaxis()->CenterTitle();
03301   hWinSigLin->SetFillColor(0);
03302   hWinSigLin->SetLineWidth(3);
03303   //hWinSigLin->SetBit(TH1::kCanRebin);
03304 
03305   TH1F* hWinAdc=new TH1F("hWinAdc","hWinAdc",3000,-1,4000);
03306   hWinAdc->SetTitle("hWinAdc");
03307   hWinAdc->GetXaxis()->SetTitle("hWinAdc");
03308   hWinAdc->GetXaxis()->CenterTitle();
03309   hWinAdc->SetFillColor(0);
03310   hWinAdc->SetLineWidth(3);
03311   //hWinAdc->SetBit(TH1::kCanRebin);
03312 
03313   TH1F* hWinPe=new TH1F("hWinPe","hWinPe",3000,-1,400);
03314   hWinPe->SetTitle("hWinPe");
03315   hWinPe->GetXaxis()->SetTitle("hWinPe");
03316   hWinPe->GetXaxis()->CenterTitle();
03317   hWinPe->SetFillColor(0);
03318   hWinPe->SetLineWidth(3);
03319   //hWinPe->SetBit(TH1::kCanRebin);
03320 
03321   TH1F* hTrkVtxPl=new TH1F("hTrkVtxPl","hTrkVtxPl",501,-1,500);
03322   hTrkVtxPl->SetTitle("hTrkVtxPl");
03323   hTrkVtxPl->GetXaxis()->SetTitle("hTrkVtxPl");
03324   hTrkVtxPl->GetXaxis()->CenterTitle();
03325   hTrkVtxPl->SetFillColor(0);
03326   hTrkVtxPl->SetLineWidth(3);
03327   //hTrkVtxPl->SetBit(TH1::kCanRebin);
03328 
03329   TH1F* hTrkEndPl=new TH1F("hTrkEndPl","hTrkEndPl",501,-1,500);
03330   hTrkEndPl->SetTitle("hTrkEndPl");
03331   hTrkEndPl->GetXaxis()->SetTitle("hTrkEndPl");
03332   hTrkEndPl->GetXaxis()->CenterTitle();
03333   hTrkEndPl->SetFillColor(0);
03334   hTrkEndPl->SetLineWidth(3);
03335   //hTrkEndPl->SetBit(TH1::kCanRebin);
03336 
03337   TH1F* hTrkLengthPl=new TH1F("hTrkLengthPl","hTrkLengthPl",501,-1,500);
03338   hTrkLengthPl->SetTitle("hTrkLengthPl");
03339   hTrkLengthPl->GetXaxis()->SetTitle("hTrkLengthPl");
03340   hTrkLengthPl->GetXaxis()->CenterTitle();
03341   hTrkLengthPl->SetFillColor(0);
03342   hTrkLengthPl->SetLineWidth(3);
03343   //hTrkLengthPl->SetBit(TH1::kCanRebin);
03344 
03345   TH1F* hWinPlFromShwEnd=new TH1F
03346     ("hWinPlFromShwEnd","hWinPlFromShwEnd",1000,-500,500);
03347   hWinPlFromShwEnd->SetTitle("hWinPlFromShwEnd");
03348   hWinPlFromShwEnd->GetXaxis()->SetTitle("hWinPlFromShwEnd");
03349   hWinPlFromShwEnd->GetXaxis()->CenterTitle();
03350   hWinPlFromShwEnd->SetFillColor(0);
03351   hWinPlFromShwEnd->SetLineWidth(3);
03352   //hWinPlFromShwEnd->SetBit(TH1::kCanRebin);
03353   
03354   TH1F* hShwVtxPl=new TH1F("hShwVtxPl","hShwVtxPl",501,-1,500);
03355   hShwVtxPl->SetTitle("hShwVtxPl");
03356   hShwVtxPl->GetXaxis()->SetTitle("hShwVtxPl");
03357   hShwVtxPl->GetXaxis()->CenterTitle();
03358   hShwVtxPl->SetFillColor(0);
03359   hShwVtxPl->SetLineWidth(3);
03360   //hShwVtxPl->SetBit(TH1::kCanRebin);
03361 
03362   TH1F* hShwEndPl=new TH1F("hShwEndPl","hShwEndPl",501,-1,500);
03363   hShwEndPl->SetTitle("hShwEndPl");
03364   hShwEndPl->GetXaxis()->SetTitle("hShwEndPl");
03365   hShwEndPl->GetXaxis()->CenterTitle();
03366   hShwEndPl->SetFillColor(0);
03367   hShwEndPl->SetLineWidth(3);
03368   //hShwEndPl->SetBit(TH1::kCanRebin);
03369 
03370   TH1F* hShwLengthPl=new TH1F("hShwLengthPl","hShwLengthPl",501,-1,500);
03371   hShwLengthPl->SetTitle("hShwLengthPl");
03372   hShwLengthPl->GetXaxis()->SetTitle("hShwLengthPl");
03373   hShwLengthPl->GetXaxis()->CenterTitle();
03374   hShwLengthPl->SetFillColor(0);
03375   hShwLengthPl->SetLineWidth(3);
03376   //hShwLengthPl->SetBit(TH1::kCanRebin);
03377 
03378   TH1F* hTrkShwVtxDiff=new TH1F("hTrkShwVtxDiff","hTrkShwVtxDiff",
03379                                 1000,-500,500);
03380   hTrkShwVtxDiff->SetTitle("hTrkShwVtxDiff");
03381   hTrkShwVtxDiff->GetXaxis()->SetTitle("hTrkShwVtxDiff");
03382   hTrkShwVtxDiff->GetXaxis()->CenterTitle();
03383   hTrkShwVtxDiff->SetFillColor(0);
03384   hTrkShwVtxDiff->SetLineWidth(3);
03385   //hTrkShwVtxDiff->SetBit(TH1::kCanRebin);
03386   
03387   TH1F* hTrkShwEndDiff=new TH1F("hTrkShwEndDiff","hTrkShwEndDiff",
03388                                 1000,-500,500);
03389   hTrkShwEndDiff->SetTitle("hTrkShwEndDiff");
03390   hTrkShwEndDiff->GetXaxis()->SetTitle("hTrkShwEndDiff");
03391   hTrkShwEndDiff->GetXaxis()->CenterTitle();
03392   hTrkShwEndDiff->SetFillColor(0);
03393   hTrkShwEndDiff->SetLineWidth(3);
03394   //hTrkShwEndDiff->SetBit(TH1::kCanRebin);
03395   
03396   TH1F* hShwVtxPl2=new TH1F("hShwVtxPl2","hShwVtxPl2",501,-1,500);
03397   hShwVtxPl2->SetTitle("hShwVtxPl2");
03398   hShwVtxPl2->GetXaxis()->SetTitle("hShwVtxPl2");
03399   hShwVtxPl2->GetXaxis()->CenterTitle();
03400   hShwVtxPl2->SetFillColor(0);
03401   hShwVtxPl2->SetLineWidth(3);
03402   //hShwVtxPl2->SetBit(TH1::kCanRebin);
03403 
03404   TH1F* hShwEndPl2=new TH1F("hShwEndPl2","hShwEndPl2",501,-1,500);
03405   hShwEndPl2->SetTitle("hShwEndPl2");
03406   hShwEndPl2->GetXaxis()->SetTitle("hShwEndPl2");
03407   hShwEndPl2->GetXaxis()->CenterTitle();
03408   hShwEndPl2->SetFillColor(0);
03409   hShwEndPl2->SetLineWidth(3);
03410   //hShwEndPl2->SetBit(TH1::kCanRebin);
03411 
03412   TH1F* hShwLengthPl2=new TH1F("hShwLengthPl2","hShwLengthPl2",501,-1,500);
03413   hShwLengthPl2->SetTitle("hShwLengthPl2");
03414   hShwLengthPl2->GetXaxis()->SetTitle("hShwLengthPl2");
03415   hShwLengthPl2->GetXaxis()->CenterTitle();
03416   hShwLengthPl2->SetFillColor(0);
03417   hShwLengthPl2->SetLineWidth(3);
03418   //hShwLengthPl2->SetBit(TH1::kCanRebin);
03419 
03420   TH1F* hTrkShwVtxDiff2=new TH1F("hTrkShwVtxDiff2","hTrkShwVtxDiff2",
03421                                 1000,-500,500);
03422   hTrkShwVtxDiff2->SetTitle("hTrkShwVtxDiff2");
03423   hTrkShwVtxDiff2->GetXaxis()->SetTitle("hTrkShwVtxDiff2");
03424   hTrkShwVtxDiff2->GetXaxis()->CenterTitle();
03425   hTrkShwVtxDiff2->SetFillColor(0);
03426   hTrkShwVtxDiff2->SetLineWidth(3);
03427   //hTrkShwVtxDiff2->SetBit(TH1::kCanRebin);
03428   
03429   TH1F* hTrkShwEndDiff2=new TH1F("hTrkShwEndDiff2","hTrkShwEndDiff2",
03430                                 1000,-500,500);
03431   hTrkShwEndDiff2->SetTitle("hTrkShwEndDiff2");
03432   hTrkShwEndDiff2->GetXaxis()->SetTitle("hTrkShwEndDiff2");
03433   hTrkShwEndDiff2->GetXaxis()->CenterTitle();
03434   hTrkShwEndDiff2->SetFillColor(0);
03435   hTrkShwEndDiff2->SetLineWidth(3);
03436   //hTrkShwEndDiff2->SetBit(TH1::kCanRebin);
03437 
03438   TH1F* hDiffYEnds=new TH1F("hDiffYEnds","hDiffYEnds",201,-4.5,4.5);
03439   hDiffYEnds->SetTitle("hDiffYEnds");
03440   hDiffYEnds->GetXaxis()->SetTitle("hDiffYEnds");
03441   hDiffYEnds->GetXaxis()->CenterTitle();
03442   hDiffYEnds->SetFillColor(0);
03443   hDiffYEnds->SetLineWidth(3);
03444   //hDiffYEnds->SetBit(TH1::kCanRebin);
03445 
03446   TH1F* hDiffYEndsSR=new TH1F("hDiffYEndsSR","hDiffYEndsSR",201,-4.5,4.5);
03447   hDiffYEndsSR->SetTitle("hDiffYEndsSR");
03448   hDiffYEndsSR->GetXaxis()->SetTitle("hDiffYEndsSR");
03449   hDiffYEndsSR->GetXaxis()->CenterTitle();
03450   hDiffYEndsSR->SetFillColor(0);
03451   hDiffYEndsSR->SetLineWidth(3);
03452   //hDiffYEndsSR->SetBit(TH1::kCanRebin);
03453 
03454   TH1F* hEvtShws=new TH1F("hEvtShws","hEvtShws",100,-50,50);
03455   hEvtShws->SetTitle("hEvtShws");
03456   hEvtShws->GetXaxis()->SetTitle("hEvtShws");
03457   hEvtShws->GetXaxis()->CenterTitle();
03458   hEvtShws->SetFillColor(0);
03459   hEvtShws->SetLineWidth(3);
03460   //hEvtShws->SetBit(TH1::kCanRebin);
03461 
03462   TProfile* pEnVsPl=new TProfile("pEnVsPl","pEnVsPl",1000,-500,500);
03463   pEnVsPl->SetTitle("Signal vs Distance from Trk Vtx");
03464   pEnVsPl->GetXaxis()->SetTitle("Distance (Planes)");
03465   pEnVsPl->GetXaxis()->CenterTitle();
03466   pEnVsPl->GetYaxis()->SetTitle("Signal");
03467   pEnVsPl->GetYaxis()->CenterTitle();
03468   pEnVsPl->SetLineColor(1);
03469   pEnVsPl->SetFillColor(0);
03470   //pEnVsPl->SetBit(TH1::kCanRebin);
03471 
03472   TProfile* pEnVsPlShw=new TProfile
03473     ("pEnVsPlShw","pEnVsPlShw",1000,-500,500);
03474   pEnVsPlShw->SetTitle("Signal vs Distance from Trk Vtx (w/ shw)");
03475   pEnVsPlShw->GetXaxis()->SetTitle("Distance (Planes)");
03476   pEnVsPlShw->GetXaxis()->CenterTitle();
03477   pEnVsPlShw->GetYaxis()->SetTitle("Signal");
03478   pEnVsPlShw->GetYaxis()->CenterTitle();
03479   pEnVsPlShw->SetLineColor(1);
03480   pEnVsPlShw->SetFillColor(0);
03481   //pEnVsPlShw->SetBit(TH1::kCanRebin);
03482 
03483   TProfile* pEnVsPlNoShw=new TProfile
03484     ("pEnVsPlNoShw","pEnVsPlNoShw",1000,-500,500);
03485   pEnVsPlNoShw->SetTitle("Signal vs Distance from Trk Vtx (no shw)");
03486   pEnVsPlNoShw->GetXaxis()->SetTitle("Distance (Planes)");
03487   pEnVsPlNoShw->GetXaxis()->CenterTitle();
03488   pEnVsPlNoShw->GetYaxis()->SetTitle("Signal");
03489   pEnVsPlNoShw->GetYaxis()->CenterTitle();
03490   pEnVsPlNoShw->SetLineColor(1);
03491   pEnVsPlNoShw->SetFillColor(0);
03492   //pEnVsPlNoShw->SetBit(TH1::kCanRebin);
03493 
03494   TProfile* pEnVsPlFromShw=new TProfile
03495     ("pEnVsPlFromShw","pEnVsPlFromShw",1000,-500,500);
03496   pEnVsPlFromShw->SetTitle("Signal vs Distance from Shw End");
03497   pEnVsPlFromShw->GetXaxis()->SetTitle("Distance (Planes)");
03498   pEnVsPlFromShw->GetXaxis()->CenterTitle();
03499   pEnVsPlFromShw->GetYaxis()->SetTitle("Signal");
03500   pEnVsPlFromShw->GetYaxis()->CenterTitle();
03501   pEnVsPlFromShw->SetLineColor(1);
03502   pEnVsPlFromShw->SetFillColor(0);
03503   //pEnVsPlFromShw->SetBit(TH1::kCanRebin);
03504 
03505   Int_t evtCounter=0;
03506   Int_t oneTrackCounter=0;
03507   Int_t oneEvtPerSlcCounter=0;
03508   Int_t longTrackCounter=0;
03509   Int_t closeTimesCounter=0;
03510   Int_t endPlaneSpectCounter=0;
03511   Int_t longTrack2Counter=0;
03512 
03513   Int_t badDirCounter=0;
03514   Int_t badDirSRCounter=0;
03515   
03516   Int_t FCCounter=0;
03517   Int_t PCCounter=0;
03518   Int_t TGCounter=0;
03519   Int_t badWindowCounter=0;
03520   Int_t badRecoCounter=0;
03521   Int_t goodWindowCounter=0;
03522   Int_t badFidCounter=0;
03523   Int_t badViewsCounter=0;
03524   Int_t badEndGapsCounter=0;
03525   Int_t badTrackEndCounter=0;
03526   Int_t bad2ndContCounter=0;
03527   Int_t badDistEndCounter=0;
03528   Int_t badShwDistCounter=0;
03529   
03530   string temps="myeventsLowEn";
03531   if (fRun>10000) temps+="MC";
03532   ofstream& eventsTxtFile=*(this->OpenTxtFile(100,"myevents"));
03533   ofstream& txtFileLowEn=*(this->OpenTxtFile(100,temps.c_str()));
03534   ofstream& txtFileNonPC=*(this->OpenTxtFile(100,"myeventsNonPC"));
03535   
03536   //reconstruction object
03537   MeuReco reco;
03538   MeuCuts cuts;
03539   MeuSummary s;
03540   
03541   //variable to turn on all the useful messages if required
03542   Msg::LogLevel_t logLevel=Msg::kDebug;
03543 
03547   
03548   this->InitialiseLoopVariables();  
03549   
03550   //for(Int_t entry=55000;entry<fEntries;entry++){
03551   for(Int_t entry=0;entry<fEntries;entry++){
03552       
03553     this->SetLoopVariables(entry);
03554     
03555     MSG("MeuAnalysis",logLevel)
03556       <<"Entry="<<entry<<endl;
03557 
03558     //make a reference instead of a pointer
03559     NtpStRecord& ntp=(*fRec);
03560     const RecCandHeader& rec=ntp.GetHeader();
03561     MAXMSG("MeuAnalysis",Msg::kInfo,1)
03562       <<"First entry: run="<<rec.GetRun()
03563       <<", snarl="<<rec.GetSnarl()<<endl;
03564     
03565     TClonesArray& evtTca=(*ntp.evt);
03566     const Int_t numEvts=evtTca.GetEntriesFast();
03567 
03568     //loop over the events
03569     for (Int_t ievt=0;ievt<numEvts;++ievt) {
03570       const NtpSREvent* pevt=
03571         dynamic_cast<NtpSREvent*>(evtTca[ievt]);
03572       const NtpSREvent& evt=(*pevt);
03573 
03574       evtCounter++;
03575 
03576       //ensure there is only 1 track in the event
03577       if (evt.ntrack!=1) continue;
03578       oneTrackCounter++;
03579         
03580       //filter out events occuring in a slice with more than 1 event
03581       if (cuts.FilterBadEvtPerSlc(ntp,evt.slc,entry)) continue;
03582       oneEvtPerSlcCounter++;
03583       
03584       TClonesArray& trkTca=(*ntp.trk);
03585       const NtpSRTrack* ptrk=
03586         dynamic_cast<NtpSRTrack*>(trkTca[evt.trk[0]]);
03587       const NtpSRTrack& trk=(*ptrk);
03588 
03589       //MeuSummary s;
03590       s.Clear();
03591       s.RFid=1.0;
03592       s.DistToEdgeFid=0.4;
03593       s.Event=entry;
03594       cuts.FillSTSumDetails(ntp,s,ievt,evt.slc);
03595       //a map to hold the info for each plane in the entry
03596       map<Int_t,MeuHitInfo> plInfo;
03597       
03598       //cut on number of planes
03599       if (sqrt(pow(1.*trk.plane.beg-trk.plane.end,2))<8) continue;
03600       longTrackCounter++;
03601       
03602       MSG("MeuAnalysis",logLevel)
03603         <<"Entry="<<entry<<", run="<<rec.GetRun()
03604         <<", snarl="<<rec.GetSnarl()
03605         <<", trk.nstrips="<<trk.nstrip<<endl;
03606       
03607       if (s.Snarl==-1){
03608         MSG("MeuAnalysis",Msg::kInfo)
03609           <<"Entry="<<entry<<", run="<<rec.GetRun()
03610           <<", snarl="<<rec.GetSnarl()
03611           <<", trk.nstrips="<<trk.nstrip<<endl;
03612         cuts.PrintNtpSt(ntp);
03613       }
03614 
03615       if (cuts.FilterBadTrkTimes(ntp,s,evt)) continue;
03616       closeTimesCounter++;
03617       
03618       cuts.ExtractPlInfo(ntp,s,plInfo,evt);
03619       
03620       reco.CalcVtx(s,plInfo);
03621       
03622       hDiffYEndsSR->Fill(plInfo[s.VtxPlane].Y-
03623                          plInfo[s.EndPlane].Y);
03624 
03625       cuts.ExtractTruthInfo(ntp,s,plInfo);
03626       Bool_t mcPosDir=s.MCVtxPlane<s.MCEndPlane;
03627       Bool_t srPosDir=trk.plane.beg<trk.plane.end;
03628       Bool_t posDir=s.VtxPlane<s.EndPlane;
03629       
03630       if (posDir!=mcPosDir && s.SimFlag==SimFlag::kMC){
03631         badDirCounter++;
03632         MAXMSG("MeuAnalysis",Msg::kVerbose,200)
03633           <<"Direction wrong! Trk: truth "
03634           <<s.MCVtxPlane<<" -> "<<s.MCEndPlane
03635           <<", reco "<<s.VtxPlane<<" -> "<<s.EndPlane
03636           //<<", y="<<s.VtxY<<" -> "<<s.EndY
03637           <<endl;
03638       }
03639       if (srPosDir!=mcPosDir && s.SimFlag==SimFlag::kMC){
03640         badDirSRCounter++;
03641         MAXMSG("MeuAnalysis",Msg::kVerbose,100)
03642           <<"Direction wrong in SR! Trk: truth "
03643           <<s.MCVtxPlane<<" -> "<<s.MCEndPlane
03644           <<", SR reco "<<trk.plane.beg<<" -> "<<trk.plane.end
03645           <<endl;
03646       }
03647 
03648       //check if track stops in spectrometer
03649       if (s.EndPlane>120) continue;
03650       endPlaneSpectCounter++;
03651       
03652       MSG("MeuAnalysis",logLevel)<<"CalcFCPC..."<<endl;
03653       reco.CalcFCPC(s,plInfo);
03654       
03655       //recalculate the vtx of the tracks
03656       //coming in from the spectrometer
03657       MSG("MeuAnalysis",logLevel)<<"CalcVtxSpecialND..."<<endl;
03658       reco.CalcVtxSpecialND(s,plInfo);
03659 
03660       //cut on number of planes with potential new vertex
03661       if (sqrt(pow(1.*s.VtxPlane-s.EndPlane,2))<8) continue;
03662       longTrack2Counter++;
03663         
03664       Bool_t TG=!s.FC && !s.PC;
03665       
03666       if (s.PC) PCCounter++;
03667       else if (s.FC) FCCounter++;
03668       else if (TG) TGCounter++;
03669       if (s.FC && s.PC) MSG("MeuAnalysis",Msg::kError)
03670         <<"Both FC and PC!"<<endl;
03671       
03672       if (s.PC) {
03673         
03674         Int_t diffAtEnd=-1;
03675         Int_t diffAtVtx=-1;
03676         if (!cuts.CheckTrkViews(s,&diffAtVtx,&diffAtEnd)) {
03677           badViewsCounter++;
03678           continue;
03679         }
03680         
03681         MSG("MeuAnalysis",logLevel)<<"Found PC event..."<<endl;
03682         Bool_t goodReco=reco.Reconstruct(s,plInfo);
03683         if (goodReco){
03684           MSG("MeuAnalysis",logLevel)<<"Recalibrating..."<<endl;
03685           this->Recalibrate(ntp,s,plInfo,evt);
03686           
03687           MSG("MeuAnalysis",logLevel)<<"CalcWindow..."<<endl;
03688           reco.CalcWindow(s,plInfo);
03689           
03690           //recalculate the containment!
03691           s.DistToEdgeFid=0.35;//allow a slight shift in position
03692           reco.CalcFCPC(s,plInfo);
03693           
03694           reco.CalcStripDists(s,plInfo);
03695           
03696           MSG("MeuAnalysis",logLevel)<<"CheckContainment..."<<endl;
03697           //fill the tree once event has been analysed
03698           if (s.WinSigMap<=0) {
03699             badWindowCounter++;
03700           }
03701           else if (!reco.CheckWinContainment(s,plInfo)){
03702             badFidCounter++;
03703             //MSG("MeuAnalysis",Msg::kInfo)<<"Outside fid:"<<endl;
03704           }
03705           else if (cuts.FilterBadTrackEnd(s,plInfo)){
03706             badTrackEndCounter++;
03707           }
03708           else if (cuts.FilterBadEndGaps(s,plInfo)){
03709             badEndGapsCounter++;
03710           }
03711           else if (!s.PC){
03712             bad2ndContCounter++;
03713             txtFileNonPC<<entry<<endl;
03714 
03715             Bool_t print=false;
03716             if (print){
03717               MAXMSG("MeuAnalysis",Msg::kDebug,50)
03718                 <<endl<<endl<<"Recalculated containment and not PC!"
03719                 <<" PC="<<s.PC<<", FC="<<s.FC
03720                 <<" entry="<<entry
03721                 <<", snarl="<<s.Snarl<<", run="<<s.Run<<endl
03722                 <<"Low En="<<s.MCLowEn<<" high="<<s.MCHighEn<<endl
03723                 <<endl;
03724               MAXMSG("MeuAnalysis",Msg::kDebug,50)
03725                 <<"vtx Pl="<<s.VtxPlane<<", end Pl="<<s.EndPlane
03726                 <<", dVtx="<<s.VtxDistToEdge
03727                 <<", dEnd="<<s.EndDistToEdge<<endl;
03728               
03729               for (map<Int_t,MeuHitInfo>::const_iterator it=
03730                      plInfo.begin();
03731                    it!=plInfo.end();++it){
03732                 const MeuHitInfo& c=it->second;
03733                 MAXMSG("MeuAnalysis",Msg::kInfo,500)
03734                   <<"pl="<<c.Plane<<", z="<<c.Z
03735                   <<", x="<<c.X<<", y="<<c.Y
03736                   <<"  t="<<c.TPos<<", l="<<c.LPos<<endl;
03737               }          
03738             }
03739           }
03740           else if (cuts.FilterBadDistEndStrip(s,plInfo)){
03741             badDistEndCounter++;
03742           }
03743           else {
03744             hTrkVtxPl->Fill(s.VtxPlane);
03745             hTrkEndPl->Fill(s.EndPlane);
03746             hTrkLengthPl->Fill(s.EndPlane-s.VtxPlane);
03747             hDiffYEnds->Fill(s.VtxY-s.EndY);
03748             hEvtShws->Fill(evt.nshower);
03749 
03750             Int_t shwVtxPl=999;
03751             Int_t shwEndPl=-999;
03752             
03753             //loop over showers in evt
03754             TClonesArray& shwTca=(*ntp.shw);
03755             for(int i=0;i<evt.nshower;i++) {    
03756               const NtpSRShower* pshw=
03757                 dynamic_cast<NtpSRShower*>(shwTca[evt.shw[i]]);
03758               const NtpSRShower& shw=(*pshw);
03759               
03760               hShwVtxPl->Fill(shw.plane.beg);
03761               hShwEndPl->Fill(shw.plane.end);
03762               hShwLengthPl->Fill(shw.plane.end-shw.plane.beg);
03763 
03764               hTrkShwVtxDiff->Fill(s.VtxPlane-shw.plane.beg);
03765               hTrkShwEndDiff->Fill(s.EndPlane-shw.plane.end);
03766 
03767               //define a vtx shower to be one that starts within 5 pl
03768               Bool_t vtxShw=abs(s.VtxPlane-shw.plane.beg)<=5;
03769               
03770               if (vtxShw){
03771                 if (shw.plane.beg<shwVtxPl) shwVtxPl=shw.plane.beg;
03772                 if (shw.plane.end>shwEndPl) shwEndPl=shw.plane.end;
03773               }
03774             }
03775 
03776             Bool_t vtxShw=false;
03777             if (shwEndPl!=-999 && shwVtxPl!=999) {
03778               vtxShw=true;
03779               hShwVtxPl2->Fill(shwVtxPl);
03780               hShwEndPl2->Fill(shwEndPl);
03781               hShwLengthPl2->Fill(shwEndPl-shwVtxPl);
03782               hTrkShwVtxDiff2->Fill(s.VtxPlane-shwVtxPl);
03783               hTrkShwEndDiff2->Fill(s.EndPlane-shwEndPl);
03784             }
03785 
03786             //fill profiles
03787             typedef map<Int_t,MeuHitInfo>::const_iterator plIt;
03788             for (plIt it=plInfo.begin();it!=plInfo.end();++it){
03789               const MeuHitInfo& cp=it->second;
03790 
03791               //check that plane is within track (assume forward going)
03792               if (cp.Plane<s.VtxPlane ||
03793                   cp.Plane>s.EndPlane) continue;
03794 
03795               Int_t plFromVtx=cp.Plane-s.VtxPlane;
03796               pEnVsPl->Fill(plFromVtx,cp.SigMap);
03797               
03798               if (vtxShw) {
03799                 pEnVsPlShw->Fill(plFromVtx,cp.SigMap);
03800                 Int_t plFromShwEnd=cp.Plane-shwEndPl;
03801                 if (plFromShwEnd>=0){
03802                   pEnVsPlFromShw->Fill(plFromShwEnd,cp.SigMap);
03803                 }
03804               }
03805               else pEnVsPlNoShw->Fill(plFromVtx,cp.SigMap);
03806             }
03807 
03808             if (vtxShw){
03809               Int_t winPlFromShwEnd=s.WinVtxSidePl-shwEndPl;
03810               hWinPlFromShwEnd->Fill(winPlFromShwEnd);
03811               
03812               //require that shw stops 10 planes before window starts
03813               if (winPlFromShwEnd<10 && s.TrigSrc>100) {
03814                 badShwDistCounter++;
03815                 continue;
03816               }
03817             }
03818             
03819 
03823             goodWindowCounter++;
03824             
03825             cuts.ExtractTruthInfo(ntp,s,plInfo);
03826             
03827             cuts.FillSTSumRecoDetails(s,plInfo);
03828             
03829             hWinSigMap->Fill(s.WinSigMap);
03830             hWinSigCor->Fill(s.WinSigCor);
03831             hWinSigLin->Fill(s.WinSigLin);
03832             hWinAdc->Fill(s.WinAdc);
03833             hWinPe->Fill(s.WinPe);
03834             
03835             eventsTxtFile<<entry<<endl;
03836             MAXMSG("MeuAnalysis",Msg::kInfo,10)
03837               <<"Good window: run="
03838               <<rec.GetRun()<<", snarl="<<rec.GetSnarl()
03839               <<endl;
03840             
03841             if (s.MCLowEn>0.5*Munits::GeV &&
03842                 s.SimFlag==SimFlag::kMC) {
03843               txtFileLowEn<<entry<<endl;
03844               MSG("MeuAnalysis",Msg::kInfo)
03845                 <<"Low En="<<s.MCLowEn<<" high="<<s.MCHighEn
03846                 <<", vtxPl="<<s.VtxPlane<<", endPl="<<s.EndPlane
03847                 <<endl<<"entry="<<entry
03848                 <<", snarl="<<s.Snarl<<", run="<<s.Run<<endl;
03849             }
03850 
03851             cuts.FillEnVsDist(s,plInfo);
03852             cuts.FillEnVsDist2(s,plInfo);
03853             cuts.FillTimeHistos(ntp,evt,s);
03854           }
03855         }
03856         else {
03857           badRecoCounter++;
03858           MSG("MeuAnalysis",logLevel)<<"Bad reco..."<<endl;
03859         }
03860       }
03861     }//end of loop over events
03862   }//end of for                                       
03863     
03867 
03868   MSG("MeuAnalysis",Msg::kInfo)<<"Finished main loop"<<endl;
03869 
03870   Double_t quantile=0.5;//quantile to compute
03871 
03872   Double_t meuWinSigMap=-1;
03873   hWinSigMap->GetQuantiles(1,&meuWinSigMap,&quantile);
03874   Double_t errWinSigMap=hWinSigMap->GetRMS()/
03875     sqrt(hWinSigMap->GetEntries());
03876   Double_t relErrWinSigMap=errWinSigMap/meuWinSigMap;
03877   
03878   Double_t meuWinSigCor=-1;
03879   hWinSigCor->GetQuantiles(1,&meuWinSigCor,&quantile);
03880   Double_t errWinSigCor=hWinSigCor->GetRMS()/
03881     sqrt(hWinSigCor->GetEntries());
03882   Double_t relErrWinSigCor=errWinSigCor/meuWinSigCor;
03883 
03884   Double_t meuWinSigLin=-1;
03885   hWinSigLin->GetQuantiles(1,&meuWinSigLin,&quantile);
03886   Double_t errWinSigLin=hWinSigLin->GetRMS()/
03887     sqrt(hWinSigLin->GetEntries());
03888   Double_t relErrWinSigLin=errWinSigLin/meuWinSigLin;
03889 
03890   Double_t meuWinAdc=-1;
03891   hWinAdc->GetQuantiles(1,&meuWinAdc,&quantile);
03892   Double_t errWinAdc=hWinAdc->GetRMS()/
03893     sqrt(hWinAdc->GetEntries());
03894   Double_t relErrWinAdc=errWinAdc/meuWinAdc;
03895 
03896   Double_t meuWinPe=-1;
03897   hWinPe->GetQuantiles(1,&meuWinPe,&quantile);
03898   Double_t errWinPe=hWinPe->GetRMS()/
03899     sqrt(hWinPe->GetEntries());
03900   Double_t relErrWinPe=errWinPe/meuWinPe;
03901 
03902   hWinPe->Fit("gaus");
03903   TF1* fPe=hWinPe->GetFunction("gaus");
03904   Double_t errWinPe2=fPe->GetParameter(2)/sqrt(hWinPe->GetEntries());
03905   
03906 //  cout<<"RESOLUTIONS (from Gaussian fit):"<<endl
03907   //    <<"ADCs:   "<<fAdc->GetParameter(2)/fAdc->GetParameter(1)<<endl
03908   //  <<"PEs:    "<<fPe->GetParameter(2)/fPe->GetParameter(1)<<endl
03909   
03910   MSG("MeuAnalysis",Msg::kInfo) 
03911     <<endl<<"*** Run Summary ***"<<endl
03912     //<<"Raw MEU="<<meu<<" +/- "<<err<<" ("<<100.*relErr<<"%)"<<endl
03913     <<"Raw MEU (WinSigMap)="<<meuWinSigMap<<" +/- "<<errWinSigMap
03914     <<" ("<<100.*relErrWinSigMap<<"%)"<<endl
03915     <<"Raw MEU (WinSigCor)="<<meuWinSigCor<<" +/- "<<errWinSigCor
03916     <<" ("<<100.*relErrWinSigCor<<"%)"<<endl
03917     <<"Raw MEU (WinSigLin)="<<meuWinSigLin<<" +/- "<<errWinSigLin
03918     <<" ("<<100.*relErrWinSigLin<<"%)"<<endl
03919     <<"Raw MEU (WinAdc)   ="<<meuWinAdc<<" +/- "<<errWinAdc
03920     <<" ("<<100.*relErrWinAdc<<"%)"<<endl
03921     <<"Raw MEU (WinPe)    ="<<meuWinPe<<" +/- "<<errWinPe
03922     <<" ("<<100.*relErrWinPe<<"%), fitErr="<<errWinPe2<<endl;
03923     //<<endl
03924     /*<<"Not LI="<<notLICounter<<endl
03925       <<"Not MM="<<notMMCounter<<endl
03926       <<"Good track="<<goodTrackCounter<<endl
03927       <<"Good length="<<goodLengthCounter<<endl
03928       <<"Good end & vtx="<<goodEndAndVtxCounter<<endl
03929       <<"Good XY reco="<<goodXYCounter<<endl
03930       <<endl
03931       <<"Good containmentPC="<<containmentCounter<<endl
03932       <<"  both SM hit (when PC)="<<bothSMCounter
03933       <<"   ("<<bothSM<<"% of PC)"<<endl
03934       <<"  coil hit (when PC)   ="<<coilHitCounter
03935       <<"   ("<<coil_PC<<"% of PC)"<<endl
03936       <<"  close LI (when PC)   ="<<closeLICounter
03937       <<"   ("<<closeLI<<"% of PC)"<<endl
03938       <<"  dead Chips (when PC) ="<<deadChipsCounter
03939       <<"   ("<<deadChips<<"% of PC)"<<endl
03940       <<"  boundary TF (when PC)="<<boundaryTFCounter
03941       <<"   ("<<boundaryTF<<"% of PC)"<<endl
03942       <<"  bad views (when PC)  ="<<badViewsCounter
03943       <<"   ("<<badViews<<"% of PC)"<<endl
03944       <<"  bad trk end (when PC)="<<badTrkEndCounter
03945       <<"   ("<<badTrkEnd<<"% of PC)"<<endl
03946       <<"  bad trk gap (when PC)="<<badTrkGapsCounter
03947       <<"   ("<<badTrkGaps<<"% of PC)"<<endl
03948       <<"  bad traceZ (when PC) ="<<badTraceZCounter
03949       <<"   ("<<badTraceZ<<"% of PC)"<<endl
03950       <<"  bad window (when PC) ="<<badWindowCounter
03951       <<"   ("<<badWin<<"% of PC)"<<endl
03952       <<"Good containmentFC="<<containmentCounterFC<<endl
03953       <<"  both SM hit (when FC)="<<bothSMCounterFC
03954       <<"   ("<<bothSMFC<<"% of FC)"<<endl
03955       <<"  coil hit (when FC)   ="<<coilHitCounterFC
03956       <<"   ("<<coil_FC<<"% of FC)"<<endl
03957       <<"Good containmentTG="<<containmentCounterTG<<endl
03958       <<"  both SM hit (when TG)="<<bothSMCounterTG
03959       <<"   ("<<bothSMTG<<"% of TG)"<<endl
03960       <<"  coil hit (when TG)   ="<<coilHitCounterTG
03961       <<"   ("<<coil_TG<<"% of TG)"<<endl
03962       <<endl
03963       <<"Bad charge pos moment="<<momentCounter<<endl
03964       <<"Too high StripsPerPlane="<<stPerPlCounter<<endl
03965       <<"  Either of above two="<<stPlOrMomentCounter<<endl
03966       <<endl
03967       <<"MC Info:"<<endl
03968       <<"  Muon lowest energy >0.6 GeV="<<lowEnCounter<<endl
03969       <<"  Muon lowest energy >5.0 GeV="<<lowEn5Counter<<endl
03970       <<endl
03971     */
03972 
03973   Int_t totalBadPCEvents=badWindowCounter+badFidCounter+
03974     badRecoCounter+badViewsCounter+
03975     badTrackEndCounter+badEndGapsCounter+bad2ndContCounter+
03976     badDistEndCounter+badShwDistCounter;
03977   
03978   MSG("MeuAnalysis",Msg::kInfo)
03979     <<"Total Entries (snarls)="<<fEntries<<endl
03980     <<"Total Evts="<<evtCounter<<endl
03981     <<"  One track      ="<<oneTrackCounter<<endl
03982     <<"  One evt in slc ="<<oneEvtPerSlcCounter<<endl
03983     <<"  Long track     ="<<longTrackCounter<<endl
03984     <<"  Close trk times="<<closeTimesCounter<<endl
03985     <<"  End plane cal  ="<<endPlaneSpectCounter<<endl
03986     <<"  Long track 2   ="<<longTrack2Counter<<endl
03987     <<"  Bad dir (MC)   ="<<badDirCounter<<endl
03988     <<"  Bad dir SR (MC)="<<badDirSRCounter<<endl
03989     <<"Containment summary:"<<endl
03990     <<"  PC="<<PCCounter<<endl
03991 //    <<"     (PC - all cuts="
03993     //closeLICounter-deadChipsCounter-boundaryTFCounter-
03994     //badViewsCounter-badTrkEndCounter-badTrkGapsCounter-badTraceZCounter-
03995     //badWindowCounter<<")"<<endl
03996     <<"  FC="<<FCCounter<<endl//<<"     (FC-coil-bothSM="
03997     //<<FCCounter-coilHitCounterFC-bothSMCounterFC<<")"<<endl
03998     <<"  TG="<<TGCounter<<endl//<<"     (TG-coil-bothSM="
03999 //    <<TGCounter-coilHitCounterTG-bothSMCounterTG<<")"<<endl
04000     <<"  FC+PC+TG="<<FCCounter+PCCounter+TGCounter<<endl
04001     <<endl
04002     <<"Total PC with good track window = "<<goodWindowCounter<<endl
04003     <<"              bad track window  = "<<badWindowCounter<<endl
04004     <<"              bad fiducial vol  = "<<badFidCounter<<endl
04005     <<"              bad reco          = "<<badRecoCounter<<endl
04006     <<"              bad views         = "<<badViewsCounter<<endl
04007     <<"              bad track end     = "<<badTrackEndCounter<<endl
04008     <<"              bad track end gaps= "<<badEndGapsCounter<<endl
04009     <<"              bad 2nd contain.  = "<<bad2ndContCounter<<endl
04010     <<"              bad dist strip end= "<<badDistEndCounter<<endl
04011     <<"              bad dist to shw   = "<<badShwDistCounter<<endl
04012     <<" Total bad                      = "<<totalBadPCEvents<<endl
04013     <<" Total bad + good               = "
04014     <<totalBadPCEvents+goodWindowCounter<<endl
04015     <<endl;
04016 
04017   if (fRec) delete fRec;
04018 
04019   MSG("MeuAnalysis",Msg::kInfo) 
04020     <<" ** Finished SpillPlots method **"<<endl;
04021 }

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

Definition at line 5018 of file MeuAnalysis.cxx.

References Msg::kInfo, MSG, MeuHistos::PrintMeuValues(), MeuCounter::PrintNtpSt(), and MeuSummaryWriter::SummaryTreeFinish().

Referenced by MeuCalModule::EndJob(), and MakeSummaryTreeWithNtpSt().

05021 {
05022   //this stores the pointers the first time the function is called
05023   static MeuCounter& cnt=*pcnt;
05024   static MeuSummaryWriter& summaryOut=*psOut;
05025   
05026   //this section runs methods on objects at the end of job
05027   if (finish){
05028     MSG("MeuAnalysis",Msg::kInfo) 
05029       <<endl<<"*** Run Summary ***"<<endl;
05030     MeuHistos histos;
05031     histos.PrintMeuValues();//this looks up TObjects by name 
05032     
05033     //print out the counter information
05034     cnt.PrintNtpSt();
05035     
05036     //finish the summary tree and output it
05037     summaryOut.SummaryTreeFinish();    
05038   }
05039 }

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, Msg::kInfo, MAXMSG, MSG, NtpSREventSummary::ndigit, NtpSRTrack::nstrip, NtpSREventSummary::ntrack, NtpSRStrip::ph1, NtpSRStrip::plane, NtpSRTrack::range, NtpSRPulseHeight::sigcor, NtpStRecord::stp, NtpSRTrack::stp, NtpSRTrack::stpz, NtpSRStrip::strip, NtpSRStrip::tpos, NtpStRecord::trk, NtpSRTrack::vtx, and NtpSRVertex::z.

00772 {
00773   MSG("MeuAnalysis",Msg::kInfo) 
00774     <<" ** Running Test method... **"<<endl;
00775 
00779   
00780   //this->InitialiseLoopVariables();  
00781   
00782   //  for(Int_t event=0;event<fEntries;event++){
00783     
00784   //this->SetLoopVariables(event);
00785 
00786   for (int i=0;i<5;i++) {
00787     fChain->GetEntry(i);
00788     NtpStRecord& ntpst=(*fRec);   // Make a reference instead of a pointer
00789     
00790     TClonesArray& tcaTk=(*ntpst.trk);
00791     Int_t numTrks=tcaTk.GetEntries();
00792    
00793     MSG("MeuAnalysis",Msg::kInfo)
00794       <<"numTrks="<<numTrks
00795       <<", ndigit="<<fRec->evthdr.ndigit
00796       <<", nTrack="<<fRec->evthdr.ntrack<<endl;
00797 
00798     //Loop over tracks
00799     for (Int_t itrk=0;itrk<numTrks;itrk++){
00800       const NtpSRTrack* ptrk=
00801         dynamic_cast<NtpSRTrack*>(tcaTk[itrk]);
00802       const NtpSRTrack& trk=(*ptrk);
00803       cout<<"range="<<trk.range<<" nstrip="<<trk.nstrip<<" ds="<<trk.ds
00804           <<" fidend.dr,dz="<<trk.fidend.dr<<" "<<trk.fidend.dz<<" end-vtx.z:"
00805           <<trk.end.z-trk.vtx.z<<" "<<trk.stpz[1]-trk.stpz[0]<<endl;
00806 
00807 
00808       TClonesArray& tcaStp=(*ntpst.stp);
00809       Int_t numStps=tcaStp.GetEntries();
00810       MSG("MeuAnalysis",Msg::kInfo)
00811         <<"numStps="<<numStps<<endl;
00812       
00813       for (Int_t i=0;i<trk.nstrip;i++){
00814 
00815         //check for bug where strip index is -1
00816         if (trk.stp[i]<0) {
00817           MAXMSG("MeuAnalysis",Msg::kInfo,500)
00818             <<"Skipping strip with trk.stp[i]="<<trk.stp[i]<<endl;
00819           continue;
00820         }
00821 
00822         //const NtpSRStrip* pstp=dynamic_cast<NtpSRStrip*>(tcaStp[i]);
00823         const NtpSRStrip* pstp=dynamic_cast<NtpSRStrip*>(tcaStp[trk.stp[i]]);
00824 
00825         const NtpSRStrip& stp=(*pstp);
00826         MSG("MeuAnalysis",Msg::kInfo)
00827           <<"Strip="<<stp.strip<<", pl="<<stp.plane
00828           <<", sigCor="<<stp.ph1.sigcor<<", tpos="<<stp.tpos<<endl;
00829       }
00830 
00831     }
00832 
00833 
00834     TClonesArray& tcaStp=(*ntpst.stp);
00835     Int_t numStps=tcaStp.GetEntries();
00836     MSG("MeuAnalysis",Msg::kInfo)
00837       <<endl<<"2nd: numStps="<<numStps<<endl;
00838     
00839     for (Int_t i=0;i<numStps;i++){      
00840       const NtpSRStrip* pstp=dynamic_cast<NtpSRStrip*>(tcaStp[i]);
00841       const NtpSRStrip& stp=(*pstp);
00842       MSG("MeuAnalysis",Msg::kInfo)
00843         <<"Strip="<<stp.strip<<", pl="<<stp.plane<<endl;
00844     }
00845     
00846     //    for (Int_t nsloop=0;nsloop<trk.nstrip;nsloop++) {
00848     //const NtpSRStrip* pstp=
00849     //dynamic_cast<NtpSRStrip*>(tcaStp[trk.stp[ns]]);
00850     //const NtpSRStrip& stp=(*pstp);
00851 
00852     //      if (stp.plane) ;   // Avoid compiler message ifndef PRINTLOTS
00853     //float du = trk.stpu[ns]-uprev;
00854     //float dv = trk.stpv[ns]-vprev;
00855     //float dz = trk.stpz[ns]-zprev;
00856     //cout<<" "<<ns<<" "<<du<<" "<<dv<<" "<<dz
00857     //    <<" uv:"<<trk.stpu[ns] <<" "<<trk.stpv[ns]
00858     //    <<" xyz:"<<trk.stpx[ns] <<" "<<trk.stpy[ns]<<" "<<trk.stpz[ns]
00859     //    <<" ds:"<<trk.stpds[ns]<<" stp: "<<trk.stp[ns]
00860     //    <<" "<<stp.plane<<endl;
00861     //  printf("%2d%8.4f%8.4f%8.4f%8.4f%2d uv:%8.4f%8.4f"
00862     //       " xyz:%8.4f%8.4f%8.4f%8.4f"
00863     //       " stp: %4d%4d\n"
00864     //       ,ns,du,dv,dz,stp.tpos,stp.planeview
00865     //       ,trk.stpu[ns],trk.stpv[ns]
00866     //       ,trk.stpx[ns],trk.stpy[ns],trk.stpz[ns]
00867     //       ,trk.stpds[ns],trk.stp[ns],stp.plane);
00868       
00869     //uprev = trk.stpu[ns];
00870     //vprev = trk.stpv[ns];
00871     //zprev = trk.stpz[ns];
00872     //}// stp loop
00873 
00874   }//end of for                                       
00875   
00879 
00880   MSG("MeuAnalysis",Msg::kInfo)<<"Finished main loop"<<endl;
00881 
00882   MSG("MeuAnalysis",Msg::kInfo) 
00883     <<" ** Finished Test method **"<<endl;
00884 }

void MeuAnalysis::TestEventLoop (  ) 

Definition at line 888 of file MeuAnalysis.cxx.

References NtpStRecord::evt, fEntries, MeuCuts::FilterBadEvtPerSlc(), InitialiseLoopVariables(), it, Msg::kDebug, Msg::kInfo, MAXMSG, MSG, NtpSREvent::nshower, NtpSREvent::nstrip, NtpSRTrack::nstrip, NtpSRShower::nstrip, NtpSRSlice::nstrip, NtpSREvent::ntrack, NtpSRStrip::ph0, NtpSRStrip::ph1, NtpSRPulseHeight::raw, SetLoopVariables(), NtpSREvent::shw, NtpStRecord::shw, NtpSREvent::slc, NtpStRecord::slc, NtpSREvent::stp, NtpStRecord::stp, NtpSRSlice::stp, NtpSRTrack::stp, NtpSRStrip::time0, NtpStRecord::trk, and NtpSREvent::trk.

00889 {
00890   MSG("MeuAnalysis",Msg::kInfo) 
00891     <<" ** Running TestEventLoop method... **"<<endl;
00892 
00893   MeuCuts cuts;
00894   
00898   
00899   this->InitialiseLoopVariables();  
00900   
00901   for(Int_t entry=0;entry<fEntries;entry++) {
00902   //for(Int_t entry=0;entry<5000;entry++) {
00903     
00904     this->SetLoopVariables(entry);
00905 
00906     NtpStRecord& ntp=(*fRec);
00907 
00908     TClonesArray& evtTca=(*ntp.evt);
00909     const Int_t numEvts=evtTca.GetEntriesFast();
00910     
00911     TClonesArray& trkTca=(*ntp.trk);
00912     const Int_t numTrks=trkTca.GetEntriesFast();
00913 
00914     TClonesArray& shwTca=(*ntp.shw);
00915     const Int_t numShws=shwTca.GetEntriesFast();
00916 
00917     TClonesArray& stpTca=(*ntp.stp);
00918     const Int_t numStps=stpTca.GetEntries();
00919 
00920     TClonesArray& slcTca=(*ntp.slc);
00921     const Int_t numSlcs=slcTca.GetEntries();
00922     
00923     Float_t totalSnarlAdc=0;
00924     Float_t totalEvtAdc=0;
00925     Float_t totalSlcAdc=0;
00926     Float_t thisSlcAdc=0;
00927     
00928     for (Int_t i=0;i<numStps;i++) {
00929       const NtpSRStrip* pstp=
00930         dynamic_cast<NtpSRStrip*>(stpTca[i]);
00931       const NtpSRStrip& stp=(*pstp);
00932       totalSnarlAdc+=stp.ph0.raw+stp.ph1.raw;
00933     }
00934     
00935     if (numEvts==0 && numTrks>0) {
00936       cout<<"Ahhh, num events in tca="<<numEvts
00937           <<", trks="<<numTrks<<endl;
00938     }
00939     
00940     if (numEvts>1) {
00941       MAXMSG("MeuAnalysis",Msg::kInfo,10)
00942         <<"High num evts="<<numEvts<<", trks="<<numTrks<<endl;
00943     }
00944 
00945     if (numSlcs!=1) {
00946       MAXMSG("MeuAnalysis",Msg::kInfo,10)
00947         <<"High num slices in tca="<<numSlcs
00948         <<", trks="<<numTrks<<endl;
00949       //continue;
00950     }
00951     
00952     MAXMSG("MeuAnalysis",Msg::kInfo,100)
00953       <<"Num events in tca="<<numEvts
00954       <<", slcs="<<numSlcs<<", trks="<<numTrks<<endl;
00955     
00956     Int_t slcStps=0;
00957     
00958     //loop over the slices
00959     for (Int_t islc=0;islc<numSlcs;++islc) {
00960       const NtpSRSlice* pslc=
00961         dynamic_cast<NtpSRSlice*>(slcTca[islc]);
00962       const NtpSRSlice& slc=(*pslc);
00963 
00964       MAXMSG("MeuAnalysis",Msg::kDebug,200)
00965         <<"  slice "<<islc+1<<" of "<<numSlcs<<endl;
00966       
00967       slcStps+=slc.nstrip;
00968 
00969       //loop over strips in slc
00970       for (Int_t i=0;i<slc.nstrip;++i) {
00971         const NtpSRStrip* pstp=
00972           dynamic_cast<NtpSRStrip*>(stpTca[slc.stp[i]]);
00973         const NtpSRStrip& stp=(*pstp);
00974         totalSlcAdc+=stp.ph0.raw+stp.ph1.raw;
00975       }
00976 
00977       //loop over tracks in slc  ---- can't do this!!!
00978       //for(int i=0;i<slc.ntrack;i++){    
00979       //const NtpSRTrack* ptrk=
00980       //  dynamic_cast<NtpSRTrack*>(trkTca[slc.trk[i]]);
00981       //const NtpSRTrack& trk=(*ptrk);
00982 
00983       //trkStps=trk.nstrip;
00984       //MAXMSG("MeuAnalysis",Msg::kInfo,300)
00985       //  <<"    track "<<i+1<<" of "<<slc.ntrack
00986       //  <<", trkStps="<<trkStps<<endl;
00987       //}
00988     }
00989 
00990     map<Int_t,Int_t> evtPerSlc;
00991     
00992     //loop over the events
00993     for (Int_t ntpevt=0;ntpevt<numEvts;++ntpevt) {
00994       const NtpSREvent* pevt=
00995         dynamic_cast<NtpSREvent*>(evtTca[ntpevt]);
00996       const NtpSREvent& evt=(*pevt);
00997       
00998       MAXMSG("MeuAnalysis",Msg::kInfo,100)
00999         <<"Entry "<<entry<<" has tracks="<<evt.ntrack
01000         <<", shws="<<numShws<<", slcs="<<numSlcs<<endl;
01001 
01002       //get the slice associated with this event
01003       const NtpSRSlice* pslc=
01004         dynamic_cast<NtpSRSlice*>(slcTca[evt.slc]);
01005       const NtpSRSlice& slc=(*pslc);
01006 
01007       evtPerSlc[evt.slc]++;
01008 
01009       //check that the filter works
01010       Int_t evts=-1;
01011       Bool_t goodEvt=!cuts.FilterBadEvtPerSlc(ntp,evt.slc,entry,&evts);
01012       if (goodEvt){
01013         MAXMSG("MeuAnalysis",Msg::kInfo,100)
01014           <<"TEST::entry="<<entry
01015           <<", evt="<<ntpevt<<" is good, evtsPerSlc="<<evts
01016           <<" (slc="<<evt.slc<<")"<<endl;
01017       }
01018       else{
01019         MAXMSG("MeuAnalysis",Msg::kInfo,100)
01020           <<"TEST::entry="<<entry
01021           <<", evt="<<ntpevt<<" is bad, evtsPerSlc="<<evts
01022           <<" (slc="<<evt.slc<<")"<<endl;
01023       }
01024       
01025       if (evt.ntrack!=1 || evt.nshower>1) continue;
01026 
01027       Int_t trkStps=0;
01028       Int_t shwStps=0;
01029       
01030       //loop over strips in slc
01031       for (Int_t i=0;i<slc.nstrip;++i) {
01032         const NtpSRStrip* pstp=
01033           dynamic_cast<NtpSRStrip*>(stpTca[slc.stp[i]]);
01034         const NtpSRStrip& stp=(*pstp);
01035         thisSlcAdc+=stp.ph0.raw+stp.ph1.raw;
01036       }
01037       
01038       //loop over tracks in evt
01039       for(int i=0;i<evt.ntrack;i++) {    
01040         const NtpSRTrack* ptrk=
01041           dynamic_cast<NtpSRTrack*>(trkTca[evt.trk[i]]);
01042         const NtpSRTrack& trk=(*ptrk);
01043 
01044         trkStps=trk.nstrip;
01045         
01046         for (Int_t i=0;i<trk.nstrip;i++) {
01047 
01048           //check for bug where strip index is -1
01049           if (trk.stp[i]<0) {
01050             MAXMSG("MeuAnalysis",Msg::kInfo,500)
01051               <<"Skipping strip with trk.stp[i]="<<trk.stp[i]<<endl;
01052             continue;
01053           }
01054 
01055           const NtpSRStrip* pstp=
01056             dynamic_cast<NtpSRStrip*>(stpTca[trk.stp[i]]);
01057           const NtpSRStrip& stp=(*pstp);
01058             MAXMSG("MeuAnalysis",Msg::kDebug,100)
01059               <<"strip time="<<stp.time0<<endl;
01060         }
01061       }
01062 
01063       //loop over showers in evt
01064       for(int i=0;i<evt.nshower;i++) {    
01065         const NtpSRShower* pshw=
01066           dynamic_cast<NtpSRShower*>(shwTca[evt.shw[i]]);
01067         const NtpSRShower& shw=(*pshw);
01068         
01069         shwStps=shw.nstrip;
01070       }
01071 
01072       //loop over strips in evt
01073       for (Int_t i=0;i<evt.nstrip;i++) {
01074         const NtpSRStrip* pstp=
01075           dynamic_cast<NtpSRStrip*>(stpTca[evt.stp[i]]);
01076         const NtpSRStrip& stp=(*pstp);
01077         totalEvtAdc+=stp.ph0.raw+stp.ph1.raw;
01078       }
01079       MAXMSG("MeuAnalysis",Msg::kInfo,100)
01080         <<"e="<<entry<<", evt="<<ntpevt
01081         <<", ADCs: snl="<<totalSnarlAdc
01082         <<", evt="<<totalEvtAdc
01083         <<", slc="<<totalSlcAdc
01084         <<", tslc="<<thisSlcAdc
01085         <<", #st="<<numStps<<","<<evt.nstrip<<","<<slcStps
01086         //<<", trk="<<trkStps<<", shw="<<shwStps
01087         <<endl;
01088     }
01089 
01090     //make this into a function
01091     //Bool_t FilterMultiEvtPerSlc(ntp,slc,entry)
01092     typedef map<Int_t,Int_t>::const_iterator slcIt;
01093     MAXMSG("MeuAnalysis",Msg::kInfo,30)
01094       <<"MAIN::Events per slice for entry="<<entry<<endl;
01095     for (slcIt it=evtPerSlc.begin();it!=evtPerSlc.end();++it) {
01096       MAXMSG("MeuAnalysis",Msg::kInfo,200)
01097         <<"  slice "<<it->first<<" has "<<it->second<<" event(s)"<<endl;
01098     }
01099     
01100   }//end of for                                       
01101     
01105   
01106   MSG("MeuAnalysis",Msg::kInfo)<<"Finished main loop"<<endl;
01107 
01108   MSG("MeuAnalysis",Msg::kInfo) 
01109     <<" ** Finished TestEventLoop method **"<<endl;
01110 }

void MeuAnalysis::UseAtNu (  )  [static]

Definition at line 215 of file MeuAnalysis.cxx.

00216 {
00217   //fUseAtNu=true;
00218 }

void MeuAnalysis::WriteOutHistos (  )  const

Definition at line 222 of file MeuAnalysis.cxx.

References fOutFile, Msg::kInfo, Msg::kWarning, and MSG.

Referenced by MeuCalModule::EndJob().

00223 {
00224   //write out the histos to the file, if it's open
00225   if (fOutFile){
00226     if (fOutFile->IsWritable()){
00227       fOutFile->cd();
00228 
00229       //create a tree and set branches      
00230       //TTree runInfo("runInfo","runInfo");
00231       //have to create local variables here
00232       //I think it's something to do with the fact that they are
00233       //read in on other trees in this class as well???
00234       //No, actually not because fSigCorsPerMip is not in a tree
00235       //don't understand this crazyness!
00236       //Int_t lSimFlag=0;
00237       //runInfo.Branch("SimFlag",&lSimFlag,"SimFlag/I");
00238 
00239 
00240       //fill and write      
00241       //runInfo.Fill();
00242       //runInfo.Write();
00243 
00244       MSG("MeuAnalysis",Msg::kInfo)
00245         <<"Writing histos to: "<<fOutFile->GetName()<<" ..."<<endl;
00246       fOutFile->Write();
00247       //fOutFile->Close();//this makes histos disappear from canvases
00248       //so do it in the destructor (need to make MeuAnalysis on heap)
00249     }
00250     else {
00251       MSG("MeuAnalysis",Msg::kWarning)
00252         <<"File not writable!"<<endl;
00253     }
00254   }
00255 }


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]
Float_t MeuAnalysis::fEntriesBD [private]

Definition at line 110 of file MeuAnalysis.h.

Referenced by SetRootFileObjects().

Definition at line 103 of file MeuAnalysis.h.

Referenced by SetRootFileObjects().

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

Definition at line 95 of file MeuAnalysis.h.

Referenced by MakeFileList().

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

Definition at line 96 of file MeuAnalysis.h.

Referenced by LoadTrees().

TFile* MeuAnalysis::fOutFile [private]

Definition at line 104 of file MeuAnalysis.h.

Referenced by BasicPlots(), N_1Plots(), SpillPlots(), WriteOutHistos(), and ~MeuAnalysis().

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

Definition at line 99 of file MeuAnalysis.h.

Referenced by Recalibrate(), and RecalibratePE().

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

Definition at line 98 of file MeuAnalysis.h.

Referenced by Recalibrate(), and RecalibrateSigCor().

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

Definition at line 97 of file MeuAnalysis.h.

Referenced by Recalibrate(), and RecalibrateSigLin().

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(), MakeSummaryTree(), and SetRootFileObjects().


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

Generated on 2 Nov 2017 for loon by  doxygen 1.6.1