MTimeEnergy.cc

Go to the documentation of this file.

#include "MTimeEnergy.hh"
#include "QCountPulsesData.hh"
#include <string>
#include <TFile.h>
#include <TTree.h>
#include <math.h>
#include <TH1D.h>
#include <TF1.h>
 
REGISTER_MODULE(MTimeEnergy)
using namespace std;
 
// constructor
MTimeEnergy::MTimeEnergy(QSequence * s) : QModule("TimeEnergy",s)
{
}
 
// destructor
MTimeEnergy::~MTimeEnergy()
{
}
 
// init method, here I set my variables to zero
void MTimeEnergy::Init() 
{
        fMultiple= new QMultiple();
        fFineTime= new QFineTime();
    string outputFilename = GetString("OutputFilename", "multiple.root");
    fOutputFile = new TFile(outputFilename.c_str(), "recreate");
    fTree = new TTree("MultipleTree", "Multiple Events");
    fTree->Branch("multiple_event_branch.", "Diana::QMultiple", &fMultiple);
        fCounter=0;
        fTotalCount=0;
        fLastTime=0;
        fSaturated=0;
        fOutput.open ("Rates.txt", ios::out | ios::app);
        fMuonFile.open("Muons.txt", ios::out);
        fDeltaT.open("DeltaT.txt", ios::out);
        fRun=0;
        for (int A1=1;A1<69;A1++)
        {
                PMCount[A1]=0.0;
                AMCount[A1]=0.0;
                MTime[A1]=0;
                AMTime[A1]=0;
                TimeHolder[A1]=0;
        }
}
 
 
 
 
        
 
 
// resetting muon scan parameters
 
void MTimeEnergy::CallReset(int run)
{
        for (int A1=1;A1<69;A1++)
        {
                if (MTime[A1]!=0)
                {
                        fOutput << A1 << "\t" << PMCount[A1] << "\t" << MTime[A1] << "\t" << AMCount[A1] << "\t" << fLastTime << "\t" << AMTime[A1] << "\n";
                }
        }               
        for (int A1=1;A1<69;A1++)
        {
                PMCount[A1]=0.0;
                AMCount[A1]=0.0;
                MTime[A1]=0;
                AMTime[A1]=0;
                TimeHolder[A1]=0;
        }
        fRun=run;
}
 
// doit method
QEvent* MTimeEnergy::Do(QEvent* ev) {
        const QRawEvent& CurrentRawEvent = ev->GetRawEvent();
 
// Adding an energy gate, perhaps odd features are caused by noise.
        double CurrentEnergy = ev->GetCalib().GetEnergy();
        
// Channel #
        unsigned long long CurrentChannel = CurrentRawEvent.GetPulse().GetChannelId();
        
// Only allow properly triggered events to be involved in the system
        int Trigger = CurrentRawEvent.GetPulse().GetIsHeater();
 
// Multiple pulse events in the same crystal are likely to be very poorly resolved in terms 
// of energy. They are possible multiple candidates, but should be ignored for now.
// Also ignoring all events without proper triggers
        int Pulses = ev->GetCountPulses().GetNumberOfPulses();
        
 
// default time for unwanted events, and so CurrentTime is set  
        unsigned long long CurrentTime = CurrentRawEvent.GetTime().GetFromStartRunNs();
 
// Use time correction routine for events which will be considered for inclusion.
        if (Pulses==1&&Trigger==0&&CurrentEnergy>=200&&CurrentChannel!=3)
        {
                CurrentTime = fFineTime->GetCrossingTime(ev);
        }
 
// Time corection due to varying position of ideal pulse triggers due to typical pulse rise times in channels
//      CurrentTime -= fFineTime->GetCorrection(CurrentChannel);
        CurrentTime += (unsigned long long) fFineTime->GetCorrection(CurrentChannel);
        
        if (fRun==0)
        {
                fRun=CurrentRawEvent.GetRun();
        }
        else if (fRun!=CurrentRawEvent.GetRun()) {CallReset(CurrentRawEvent.GetRun());}
        
        float Allowed = 20e6;
        
 
        
// Time correction due to changing trigger time with increasing energy
// Should no longer be necessary with better time correction from cfd
//      CurrentTime += 4.096*1e9*(ev->GetOB().GetStartTimeOffset())/512;
        
 
 
// Try to use the 'average' time of a multiple rather than the earliest detected time
        fLastTime = (unsigned long long) ( (CurrentTime + fLastTime*1.0*fMultiple->GetMultiplicity())/(1+1.0*fMultiple->GetMultiplicity()) );
        
// Check that new average time is consistent with other events
        int TimeCheck = fMultiple->Consistancy(fLastTime,Allowed);
        unsigned long long DeltaT = fLastTime-CurrentTime;
        
        
// Demand a good pulse  
        if (Pulses==1&&Trigger==0&&CurrentEnergy>=200&&CurrentChannel!=3)
        {
                if (ev->GetPreProcess().GetIsSaturatedHigh()==1) {fSaturated++;}
                if (fCounter==50)
                {
                        fDeltaT << CurrentTime-TimeHolder[CurrentChannel] << "\n";
                }
                TimeHolder[CurrentChannel]=CurrentTime;
 
                DeltaT*=DeltaT;
                DeltaT= (unsigned long long) pow(DeltaT,0.5);
// no data in active multiple (first event), no check, add current event
                if (fMultiple->GetMultiplicity() == 0)
                {
                        fLastTime=CurrentTime;
                        fMultiple->AddEvent(ev,CurrentTime);
                        fTotalCount++;
                }
// multiple member check, if succuss, add current event to new multiple
// also recording the total time splitting in the event with AccumulateTime()
// although typical time splitting for simultaneous low energy events has been 
// reduced to 2 ms, at higher energies the pulse shapes and pulse fitting is 
// not so well known, and thus time finding is less accurate. This also checks
// that the current time is acceptable to the other members of the multiple
                else if (DeltaT<Allowed&&DeltaT>-Allowed&&TimeCheck==1)
                {
                        fMultiple->AddEvent(ev,CurrentTime);
                        fTotalCount++;
                        fMultiple->AccumulateTime(DeltaT);
                }
// failed check, end current multiple, start new multiple with current event
// check current multiple for muonlike status 
                else
                {
// checks for total energy, and presence of saturated peaks
                        if (fMultiple->GetTotalEnergy()>30000&&fSaturated>0)
                        {
                                fMultiple->SetMuon();
                                for (int A2=0;A2<fMultiple->GetMultiplicity();A2++)
                                {
                                        int TaggedChannel = fMultiple->GetEventChannelID(A2);
                                        // only count earliest tag, avoids confusion over count.
                                        if (MTime[TaggedChannel]==0) {MTime[TaggedChannel]=CurrentTime;}
                                }
                                const QBaseRawEventR* HolderRawEvent = fMultiple->GetEvent(0).GetBaseRawEvent();
                                fMuonFile << MTime[fMultiple->GetEventChannelID(0)] << "\t" << HolderRawEvent->GetRun() << "\t" << fMultiple->GetMultiplicity() << "\t";
                                for (int A2=0;A2<fMultiple->GetMultiplicity();A2++)
                                {
                                        fMuonFile << fMultiple->GetEventChannelID(A2) << "\t";
                                }
                                fMuonFile << "\n";
                        }
                        
// checks for presence of events from poor channels/runs/times in the multiple and flags the multiple if they are present.
                        fMultiple->CheckFailures();
 
// erasing current multiple and starting a new one
                        fMultiple->EliminateSize();
                        fTree->Fill();
                        delete fMultiple;
                        fMultiple = new QMultiple();
                        fMultiple->AddEvent(ev,CurrentTime);
                        fTotalCount++;
                        fLastTime=CurrentTime;
                        if (MTime[CurrentChannel]==0)
                        {
                                PMCount[CurrentChannel]++;
                        }
                        else if (CurrentTime-MTime[CurrentChannel]<1e12) 
                        {
                                AMCount[CurrentChannel]++;
                                if (AMCount[CurrentChannel]==1)
                                {
                                        AMTime[CurrentChannel]=CurrentTime-MTime[CurrentChannel];
                                }
                        }
                        fSaturated=0;
                        fCounter++;
                        if (fCounter>51)
                        {
                                cout << ev->GetReadNumber() << "\t" << fTotalCount << "\n";
                                fCounter =0;
                        }
                }
        }
        return ev;
}
 
// done method
void MTimeEnergy::Done()
{
        for (int A1=1;A1<69;A1++)
        {
                if (MTime[A1]!=0)
                {
                        fOutput << A1 << "\t" << PMCount[A1] << "\t" << MTime[A1] << "\t" << AMCount[A1] << "\t" << fLastTime << "\t" << AMTime[A1] << "\n";
                }
        }                       
        fOutputFile = fTree->GetCurrentFile();
    fOutputFile->Write();
    fOutputFile->Close();
}