/home/alex/Work/Research/Km3NeT/PPM-DU/ppm-du-v2.3/analyser.cc

Go to the documentation of this file.

#include <fstream>
#include <vector>
#include <string>
#include <ctime>
#include <dirent.h>
#include <boost/lexical_cast.hpp>
#include <math.h>
//#include <boost/filesystem/operations.hpp>

#include <TSystem.h>
#include <TFile.h>
#include <TTree.h>
#include <TBranch.h>
#include <TLeaf.h>

#include "rawReader.h"
#include "event.h"
#include "dom.h"
#include "pmt.h"
#include "signal.h"
#include "signals.h"
#include "triggeringDOM.h"
#include "triggeredSignals.h"
#include "triggeringEvent.h"
#include "histo.h"
#include "constants.h"


using namespace std;

int
main(int argc, char* argv[])
{
  time_t time0, timeA, timeB;
  time(&time0);
  time(&timeA);
  vector<dom> vSlice;
  vector<triggeredSignals> vEvents;
  vector<int> vDomId;
  vDomId.push_back(101);
  vDomId.push_back(102);
  vDomId.push_back(103);
  vEvents.clear();
  int eventTotal = 0;
  int event2Doms = 0;
  int eventTrigger = 0;
  //const bool individualPmts = false;
  //const int refId = -1;
  TFile hfile(argv[3], "RECREATE");
  histo hResults(vDomId[domId], domId, pmtId);
  hResults.InitEventPerSlice(individualPmts, vDomId);
  hResults.InitSignal();
  hResults.InitFrequency(vDomId);
  hResults.InitDurations(vDomId);
  hResults.InitDeltaT(vDomId);
  hResults.InitDomDeltaT(vDomId);
  hResults.InitSliceTimeSeries(vDomId);
  hResults.InitCoincidenceRate(vDomId);
  long long int startRunTime = 0;
  bool startSet = false;
  bool toWrite = false;
  long long int endRunTime = 0;
  if (argc < 4) { 
    cerr << "give the id of the first and last";
    cerr << " runs you want to analyse and the root file name!" << endl;
    return 0;
  }
  int sliceTot = 0;
  int nbSlicePmt[domNb] = {0, 0, 0};
  int badSlice[domNb] = {0, 0, 0};
  for (int r = boost::lexical_cast<int>(argv[1]); r <= boost::lexical_cast<int>(argv[2]); ++r) {
    struct dirent* dirStruc;
    const string dataPath = "/home/alex/Data/Km3NeT/PPM-DU/raw/";
    hResults.SetRunId(r);
    DIR* directory = opendir(dataPath.c_str());
    string filename;
    char filenameStart[1024];
    sprintf(filenameStart, "RUN-PPM_DU-%0*d", 5, r);
    char filenameEnd[1024];
    sprintf(filenameEnd, ".rdat");
    while ((dirStruc = readdir(directory))) {
      const string dirname = string(dirStruc->d_name);
      if ((dirname.find(filenameStart) != string::npos)
          && (dirname.find(filenameEnd) != string::npos)) {
        filename = dataPath + dirname;
      }
    }
    bool goodRun = true;
    for (int g = 4; g <= argc - 1; ++g) {
      const int wrongId = boost::lexical_cast<int>(argv[g]);
      if (r == wrongId) goodRun = false;
    }
    ifstream file;
    file.open(filename.c_str(), ios::binary);
    int twoGbSteps = 1;
    if ((!file) || (goodRun == false)) {
      cerr << r << " is not the correct file name!" << endl;
      continue;
    } else {
      long long int fileLength = GetFileLength(file);
      twoGbSteps = int(double(fileLength)/2.0e9) + 1;
      if (twoGbSteps > 1) {
        cerr << "file " << filename.c_str() << " is larger than 2.0GB - "
             << twoGbSteps << " loops will be made" << endl;
      }
      long long int bitPosition = file.tellg();
      time(&timeA);
      cerr << "********** looking for the cut values **********\n";   
      for (int i = 0; i < twoGbSteps; ++i) {
        vSlice.clear();
        vSlice.resize(0);
        bitPosition = FillRaw(file, vSlice, bitPosition, toWrite);
        for (unsigned int s = 0; s < vSlice.size(); ++s) {
          const int curId = vSlice[s].GetDomHeader().GetDomId();
          if (curId == vDomId[domId]) {
            const long long int curtime = vSlice[s].GetDomHeader().GetDomTime();
            if (vSlice[s].HasPmtData(pmtId) == true) {
              pmt curPmt = vSlice[s].GetPmtData(pmtId);
              if (curPmt.HasSignal() == true) {
                const vector<signal> vSignal = curPmt.GetSignal();
                hResults.FillSignal(vSignal, curtime);
              }
            }
          }
          bool countedDom = false;
          for (unsigned int j = 0; j < vDomId.size(); ++j) {
            if (curId == vDomId[j]) {
              countedDom = true;
              break;
            }
          }
          if (vDomId.size() == 0) countedDom = false;
          if (countedDom == false) vDomId.push_back(curId);
        }
        hResults.FillSliceTimeSeries(vSlice, vDomId);
        hResults.FillEventPerSlice(vSlice, individualPmts, vDomId);
      }
    }
    file.close();
    double cutValues[domNb][pmtNb];
    hResults.DrawEventPerSlice(individualPmts, cutValues, vDomId);
    hResults.DrawSignal();
    file.open(filename.c_str(), ios::binary);
    long long int bitPosition = file.tellg();
    time(&timeB);
    cerr << "cut values are (in kHz):\n";
    for (int dom = 0; dom < domNb; ++dom) {
      for (int pmt = 0; pmt < pmtNb; ++pmt) {
        //cutValues[dom][pmt] = DBL_MAX; // to remove the cuts!
        cerr << "dom" << dom << " - pmt" << pmt << " - " << cutValues[dom][pmt]/conversion << endl;
      }   
    } 
    cerr << "time of analysis: " << difftime(timeB, timeA) << 's' << endl;
    time(&timeA);
    toWrite = false;
    vector<triggeredSignals> v2GigaResults;
    for (int i = 0; i < twoGbSteps; ++i) {
      //v2GigaResults.clear();
      cerr << "********** Filling the DOM with data - loop " << i << " **********\n";   
      vSlice.clear();
      bitPosition = FillRaw(file, vSlice, bitPosition, toWrite);
      time(&timeB);
      cerr << "time of analysis: " << difftime(timeB, timeA) << 's' << endl;
      cerr << "********** Filling histograms **********\n"; 
      //const double cutValue = cutValues[cutValues.size() - 1];
      hResults.FillFrequency(vSlice, cutValues, vDomId);
      hResults.FillDurations(vSlice, cutValues, vDomId);
      time(&timeA);
      cerr << "time of analysis: " << difftime(timeA, timeB) << 's' << endl;     
      cerr << "********** Filling events and first level trigger **********\n";
      vector<triggeredSignals> vTrig2pmtResult;
      vTrig2pmtResult.clear();
      vector<triggeredSignals> vTrigXpmtResult;
      vTrigXpmtResult.clear();
      vector<signals> vSignals;
      vSignals.clear();
      sliceTot += vSlice.size();
      for (unsigned int e = 0; e < vSlice.size(); ++e) {
        //cerr << 100.*double(e)/double(vSlice.size()) << "%        \r";
        const domHeader& domHead = vSlice[e].GetDomHeader();
        const long long int domTime = domHead.GetDomTime();
        if (startSet == false) {
          startRunTime = domTime;
          startSet = true;
        }
        if (e == vSlice.size() - 1) endRunTime = domTime;
        const int domId = domHead.GetDomId();
        int domIndex = hResults.FindDomIndex(domId, vDomId);
        hResults.SetDomId(domIndex);
        //hResults.InitSignals();
        const vector<pmt>& pmtSample = vSlice[e].GetPmtsData();
        if (!pmtSample.size()) continue;
        int sigInSlice = 0;
        for (unsigned int j = 0; j < pmtSample.size(); ++j) {
          vector<signal> sigSample = pmtSample[j].GetSignal();
          const int sigSampleSize = sigSample.size();
          sigInSlice += sigSampleSize;
          if (!sigSampleSize) continue;
          //const double cutValue = cutValues[cutValues.size() - 1];
          const int pmtId = pmtSample[j].GetPmtId();
          if (IsNoisySlice(sigSampleSize, cutValues[domIndex][pmtId])) {
            ++badSlice[domIndex];
          } else { 
            ++nbSlicePmt[domIndex];
            for (int k = 0; k < sigSampleSize; ++k) {
              const long long int start = sigSample[k].GetStart() + domTime;
              const signals oneSig(domIndex, pmtId, start, start + sigSample[k].GetDuration());
              FillSignals(oneSig, vSignals);
            }
          }
        }
        eventTotal += sigInSlice;
      }
      SortSignals(vSignals);
      //for (unsigned int sig = 0; sig < vSignals.size(); ++sig) {
      //cout << vSignals[sig].GetDomId() << ' ' << vSignals[sig].GetPmtId() << ' ' 
      //     << vSignals[sig].GetStart() << ' ' << vSignals[sig].GetStop() << endl;             
      //}
      time(&timeB);
      DoT2Pmts(vSignals, vTrig2pmtResult, coincWindowPmts);
      DoTmultiPmts(vSignals, vTrigXpmtResult, coincWindowDoms);
      int localmulti = 4;
      char f4foldname[1024];
      sprintf(f4foldname, "/home/alex/Data/Km3NeT/PPM-DU/stats/RUN-PPM_DU-%0*d-%dfold-%dns.dat",
              5, r, localmulti, coincWindowPmts);
      ofstream out4fold(f4foldname, ios::app);
      for (unsigned int i = 0; i < vTrigXpmtResult.size(); ++i) {
        const std::vector<signals>& vSignals4fold = vTrigXpmtResult[i].GetSignals();
        if (vSignals4fold.size() == unsigned(localmulti)) {
          out4fold << vSignals4fold[0].GetDomId() << ' ';
          for (unsigned int j = 0; j < vSignals4fold.size(); ++j) {
            out4fold << vSignals4fold[j].GetPmtId() << ' ' << vSignals4fold[j].GetStart()
                     << ' ' << vSignals4fold[j].GetStop() << ' ';
          }
          out4fold << endl;
        }
      }
      out4fold.close();
      //for (unsigned int i = 0; i < vTrig2pmtResult.size(); ++i) {
      //const std::vector<signals>& vSignals = vTrig2pmtResult[i].GetSignals();
      //const int pmtId = vSignals[0].GetPmtId();
      //if (vSignals.size() > 2)
      //  cout << pmtId << ' ' << vSignals.size() << endl;
      //}
      hResults.FillCoincidenceRate(vTrigXpmtResult);
      hResults.FillDeltaT(vTrig2pmtResult);
      vector<triggeredSignals> vTrig2Result;
      SelectingMultiplicity(vTrigXpmtResult, vTrig2Result, multiplicityForMuon);
      AddDom2Event(vTrig2Result, v2GigaResults);        
      //}
      hResults.DrawPmtStats(vSlice, vDomId);
      time(&timeB);
      cerr << "time of analysis: " << difftime(timeB, timeA) << 's' << endl;    
    }
    cerr << "*************** Second level trigger ****************" << endl;
    Do2DomTrigger(v2GigaResults, vEvents, domTimeWindow);
    //DoEventTrigger(v2GigaResults, vEvents, domTimeWindow);
    eventTrigger = vEvents.size();
    event2Doms = v2GigaResults.size();
    hResults.FillDomDeltaT(vEvents, vDomId);
    hResults.DrawDurations();
    hResults.DrawFrequency();
    hResults.DrawDeltaT();
    hResults.DrawSliceTimeSeries();
    hResults.DrawDomDeltaT();
    //if (eventTrigger != 0) hResults.FillDomDeltaT(vEvents, vDomId);
    cout << r << ' ' << sliceTot << " 0 " << nbSlicePmt[0] << ' ' << badSlice[0]
         << " 1 " << nbSlicePmt[1] << ' ' << badSlice[1]
         << " 2 " << nbSlicePmt[2] << ' ' << badSlice[2] << endl;
  }
  hResults.DrawCoincidenceRate();

  hfile.Close();
  cerr << "number of not noisy pmts - number of noisy pmts\n"
       << "dom 101: " << nbSlicePmt[0] << ' ' << badSlice[0] << "\n" 
       << "dom 102: " << nbSlicePmt[1] << ' ' << badSlice[1] << "\n"
       << "dom 103: " << nbSlicePmt[2] << ' ' << badSlice[2] << "\n"
       << "Total number of signals: " << eventTotal
       << " - number of events for a multiplicity of " << multiplicityForMuon << ": " << event2Doms
       << " - number of 2 Doms signals: " << eventTrigger << endl;
  time(&timeA);
  const double durationInS = double(endRunTime - startRunTime)*1e-9 + 0.134217728;
  const int runHour = int(durationInS/3600.);
  const int runMin = int((durationInS - double(3600*runHour))/60.);
  const double runSec = durationInS - double(3600*runHour) - double(60*runMin);
  cerr << "time of analysis: " << difftime(timeA, time0) << "s\n"
       << "run duration: " << runHour << "h " << runMin << "min " << runSec << "s"<< endl;
  return 0;
}