beamtime/cern-oct11/go4/RICH/TRICHProc.cxx (r4864/r3566)

Go to the documentation of this file.

/* Generated by Together */

#include "TRICHProc.h"
#include "TCBMBeamtimeEvent.h"

#include "TGo4Version.h"
#if __GO4BUILDVERSION__ > 40502
#include "go4iostream.h"
#else
#include "Riostream.h"
#endif
#include "TSystem.h"
#include "TGo4UserException.h"
#include "TGo4MbsEvent.h"
#include "roc/Message.h"
#include "roc/Board.h"

#include "TGo4Analysis.h"

#include "CbmRichRingFitterCOP.h"
#include "CbmRichRingFitterEllipseTau.h"
#include "CbmRichRingFinderHoughImpl.h"

#include "TEllipse.h"
#include "TLatex.h"
#include <sstream>
#include <fstream>

TRICHProc::TRICHProc() :
   TCBMBeamtimeProc(),
   fRocInputEvent(0),
   fOutputEvent(0),
   fHodo1(0),
   fTriglogEvent(0),
   fBeamEvent(0),
   fEventNr(0),
   fNrSED(0)
{
}


TRICHProc::TRICHProc(const char* name) :
   TCBMBeamtimeProc(name),
        fRocInputEvent(0),
        fOutputEvent(0),
    fHodo1(0),
    fTriglogEvent(0),
        fBeamEvent(0),
        fEventNr(0),
        fNrSED(0)
{
   cout << "**** TRICHProc: Create instance " << name << endl;
   
   
   fPar = (TRICHParam*) MakeParameter("RICHPar", "TRICHParam");
   
   
   if (!ReadMapmtGeometry("mapmt-geometry.txt")) {
      printf("First step of reading MAPMT geometry failed \n");
      if (!ReadMapmtGeometry("RICH/mapmt-geometry.txt")) {
         printf("second step of reading MAPMT geometry failed \n");
      }
   }
   
   hEntries2dIntegral = MakeTH2('D', "MAPMT/Entries2dIntegral", "Integral of ADC values per MAPMT channels", 32, 0., 32., 32, 0., 32., "x-pos (pixel)", "y-pos (pixel)");
   hEntries2dMean = MakeTH2('D', "MAPMT/Entries2dMean", "Mean ADC values per MAPMT channels", 32, 10., 32., 32, 0., 32., "x-pos (pixel)", "y-pos (pixel)");
   
   hEntriesperPMT = MakeTH2('I', "MAPMT/EntriesperPMT", "Number of hits per MAPMT", 4, 0., 32., 4, 0., 32., "x-pos (pixel)", "y-pos (pixel)");
   
   hADCall = MakeTH1('I', "MAPMT/ADCall", "Distribution of corrected ADC values", 100, 0., 4096.);
   hHittime = MakeTH1('I', "MAPMT/Hittime", "Distribution of difference between AUX2 and MAPMT hit", 1024, -1500, 4096.);

   for (Int_t iPMT=1; iPMT<=16; iPMT++) {
      hHitMultperPMT[iPMT]=MakeTH1('I', Form("MAPMT/HitMultperPMT%i",iPMT), Form("Hit Multiplicity on PMT %i",iPMT),50,-0.5,49.5, "x (mm)", "y (mm)");
   }
   


   // window scan of the hodoscope   
   fNofHodWindowsX = 4;
   fNofHodWindowsY = 4;
   hHodXYscan=MakeTH2('F', "MAPMT/HodXYscan", "HodXYscan", 16*fNofHodWindowsX, 0., 16.*fNofHodWindowsX, 16.*fNofHodWindowsY, 0., 16.*fNofHodWindowsY, "x-pos (pixel)", "y-pos (pixel)");
      
   hNofHitsInHodo = MakeTH1('F', "MAPMT/hNofHitsInHodo", "Number of hits in hodoscope", 4, 0, 4, "Number of hits in hodoscope");
   hHodoADCvsNofHitInEv = MakeTH2('F', "MAPMT/hHodoADCvsNofHitInEv", "ADC value in hodoscope vs. number of hits in event", 6000, 0., 6000., 60, 0., 60., "ADC in hodoscope", "Number of hits in event");
   hNotAssignedHitsXY = MakeTH2('F', "MAPMT/hNotAssignedHitsXY", "Not assigned hits", 32, 0., 208., 32, 0., 208., "x-pos (mm)", "y-pos (mm)");

   // book additional 2d histograms for position plots with condition on fiber hod
   hMAPMT_hcor2d_s1 = MakeTH2('I', "MAPMT/mapmt_detector_h1", "Number of hits per MAPMT channels, hcut1", 32, 0., 32., 32, 0., 32., "x-pos (pixel)", "y-pos (pixel)");
   hMAPMT_hcor2d_s2 = MakeTH2('I', "MAPMT/mapmt_detector_h2", "Number of hits per MAPMT channels, hcut2", 32, 0., 32., 32, 0., 32., "x-pos (pixel)", "y-pos (pixel)");
   hMAPMT_hcor2d_s3 = MakeTH2('I', "MAPMT/mapmt_detector_h3", "Number of hits per MAPMT channels, hcut3", 32, 0., 32., 32, 0., 32., "x-pos (pixel)", "y-pos (pixel)");
   hMAPMT_hcor2d_s4 = MakeTH2('I', "MAPMT/mapmt_detector_h4", "Number of hits per MAPMT channels, hcut4", 32, 0., 32., 32, 0., 32., "x-pos (pixel)", "y-pos (pixel)");
   hMAPMT_hcor2d_s5 = MakeTH2('I', "MAPMT/mapmt_detector_h5", "Number of hits per MAPMT channels, hcut5", 32, 0., 32., 32, 0., 32., "x-pos (pixel)", "y-pos (pixel)");
   
   hNofFoundRings= MakeTH1('F', "MAPMT/hNofFoundRings", "Number of found rings", 3, 0.0, 3., "Number of found rings");

   hTimespread=MakeTH1('I',"MAPMT/Timespread","Timespread of hits",50,-25,25,"time (ns)","hits");
   
   //[0] is electron, [1] is pion, [2] is muon, [3] is all
   string p[] = {"e", "pi", "mu", "all"};
   for (Int_t i = 0; i < 4; i++){
      //book circle fitting histograms
      hCircleFitRadius[i] = MakeTH1('F', (string("MAPMT/fitting/")+p[i]+string("/CircleFitRadius_")+p[i]).c_str(), "Radius of fitted rings", 501, -0.5, 100.5, "ring radius [mm]");
      hCircleFitCenter[i] = MakeTH2('F', (string("MAPMT/fitting/")+p[i]+string("/CircleFitCenter_")+p[i]).c_str(), "Center of fitted rings", 208, 0., 208., 208, 0., 208., "x-pos [mm]", "y-pos [mm]");
      hCircleFitChi2[i] = MakeTH1('F', (string("MAPMT/fitting/")+p[i]+string("/CircleFitChi2_")+p[i]).c_str(), "Chi2 of fitted rings", 200, 0., 40., "chi2");
      hCircleFitDR[i] = MakeTH1('F', (string("MAPMT/fitting/")+p[i]+string("/CircleFitDR_")+p[i]).c_str(), "dR of fitted rings", 200, -20., 20., "dR [mm]");
      
      //book ellipse fitting histograms
      hEllipseFitAaxis[i] = MakeTH1('F', (string("MAPMT/fitting/")+p[i]+string("/EllipseFitAaxis_")+p[i]).c_str(), "A axis of fitted ellipse", 501, -0.5, 100.5, "A axis [mm]");
      hEllipseFitBaxis[i] = MakeTH1('F', (string("MAPMT/fitting/")+p[i]+string("/EllipseFitBaxis_")+p[i]).c_str(), "B axis of fitted ellipse", 501, -0.5, 100.5, "B axis [mm]");
      hEllipseFitPhi[i] = MakeTH1('F', (string("MAPMT/fitting/")+p[i]+string("/EllipseFitPhi_")+p[i]).c_str(), "Phi of fitted ellipses", 400, -2., 2., "Phi [rad]");
      hEllipseFitBoverA[i] = MakeTH1('F', (string("MAPMT/fitting/")+p[i]+string("/EllipseFitBoverA_")+p[i]).c_str(), "B/A of fitted ellipses", 100, 0., 1., "B/A");
      hEllipseFitCenter[i] = MakeTH2('F', (string("MAPMT/fitting/")+p[i]+string("/EllipseFitCenter_")+p[i]).c_str(), "Center of fitted ellipses", 208, 0., 208., 208, 0., 208., "x-pos [mm]", "y-pos [mm]");
      hEllipseFitChi2[i] = MakeTH1('F', (string("MAPMT/fitting/")+p[i]+string("/EllipseFitChi2_")+p[i]).c_str(), "Chi2 of fitted ellipses", 200, 0., 40., "chi2");

      //number of hits per event for all events, for events with good circle fit, for event with good ellipse fit
      hNofHitsEvAll[i] = MakeTH1('F', (string("MAPMT/efficiency/")+p[i]+string("/NofHitsEvAll_")+p[i]).c_str(), "Nof hits in event", 60, 0, 60, "Nof hits in event");
      hNofHitsEvRingFinder[i] = MakeTH1('F', (string("MAPMT/efficiency/")+p[i]+string("/NofHitsEvRingFinder_")+p[i]).c_str(), "Nof hits in event", 60, 0, 60, "Nof hits in event");
      hNofHitsEvCircleFit[i] = MakeTH1('F', (string("MAPMT/efficiency/")+p[i]+string("/NofHitsEvCircleFit_")+p[i]).c_str(), "Nof hits in event", 60, 0, 60, "Nof hits in event");
      hNofHitsEvEllipseFit[i] = MakeTH1('F', (string("MAPMT/efficiency/")+p[i]+string("/NofHitsEvEllipseFit_")+p[i]).c_str(), "Nof hits in event", 60, 0, 60, "Nof hits in event");
      hNofHitsEvRingFinderEff[i] = MakeTH1('F', (string("MAPMT/efficiency/")+p[i]+string("/NofHitsEvRingFinderEff_")+p[i]).c_str(), "Nof hits in event", 60, 0, 60, "Nof hits in event");
      hNofHitsEvCircleFitEff[i] = MakeTH1('F', (string("MAPMT/efficiency/")+p[i]+string("/NofHitsEvCircleFitEff_")+p[i]).c_str(), "Nof hits in event", 60, 0, 60, "Nof hits in event");
      hNofHitsEvEllipseFitEff[i] = MakeTH1('F', (string("MAPMT/efficiency/")+p[i]+string("/NofHitsEvEllipseFitEff_")+p[i]).c_str(), "Nof hits in event", 60, 0, 60, "Nof hits in event");

      hNofHitsRingFinder[i] = MakeTH1('F', (string("MAPMT/efficiency/")+p[i]+string("/NofHitsRingFinder_")+p[i]).c_str(), "Nof hits in ring", 60, 0, 60, "Nof hits in ring");

      hEntries2d[i]=MakeTH2('I', (string("MAPMT/Entries2d_")+p[i]).c_str(), "Hits per Pixel", 32, 0., 32., 32, 0., 32., "x-pos [pixel]", "y-pos [pixel]");
   }

   SetParticleIdentificationConditions();
   
   hMAPMT_HitMultvsO2 = MakeTH2('I', "MAPMT/mapmt_hitmultvsO2","HitMultiplicity vs O2",100,0,400,50,-0.5,49.5,"O2 level","hitmult");
   
   // book conditions
   fMAPMT_cond = MakeWinCond("hADCall", 20., 2500., "hADCall"); //0., 1500.
   fMAPMT_cond->SetHistogram("hADCall");
   fMAPMT_tcond = MakeWinCond("hHittime", -200., -100., "hHittime"); //0., 900.
   fMAPMT_tcond->SetHistogram("hHittime");

   Int_t nofPics = 30;
   int maxNofRings = 2;
   fCirclePicture.resize(nofPics);
   fCircle.resize(nofPics);
   fCircleText.resize(nofPics);
   fEllipse.resize(nofPics);
   fEllipseText.resize(nofPics);
   hSED.resize(nofPics);
   fRingHits.resize(nofPics);
   fEventHits.resize(nofPics);
   for (Int_t i = 0; i < nofPics; i++){

      hSED[i] = MakeTH2('I', Form("MAPMT/SED/SED%i",i),Form("Single Event %i",i), 32, 0., 208., 32, 0., 208., "x-pos (mm)", "y-pos (mm)");

      fCirclePicture[i] = new TGo4Picture(Form("Event%i",i), "Circle picture");

      fCircle[i].resize(maxNofRings);
      fEllipse[i].resize(maxNofRings);
      for (int iR = 0; iR < maxNofRings; iR++){
                  fCircle[i][iR] = new TEllipse(16., 16., 10.);
                  fCircle[i][iR]->SetFillStyle(0);
                  fCircle[i][iR]->SetLineWidth(2);

                  fEllipse[i][iR] = new TEllipse(16., 16., 10.);
                  fEllipse[i][iR]->SetFillStyle(0);
                  fEllipse[i][iR]->SetLineWidth(2);
                  fEllipse[i][iR]->SetLineColor(kBlue);

              fCirclePicture[i]->AddSpecialObject(fCircle[i][iR]);
              fCirclePicture[i]->AddSpecialObject(fEllipse[i][iR]);
      }

      fCircleText[i] = new TLatex(5, 30, "(x,y,r)");
      fEllipseText[i] = new TLatex(5, 30, "(x,y,r)");

      fCirclePicture[i]->AddObject(hSED[i]);
      fCirclePicture[i]->AddSpecialObject(fCircleText[i]);
      //fCirclePicture[i]->AddSpecialObject(fCircle[i]);
      fCirclePicture[i]->AddSpecialObject(fEllipseText[i]);
      //fCirclePicture[i]->AddSpecialObject(fEllipse[i]);

      fRingHits[i].resize(fMaxNofHitsInEventDraw);
      fEventHits[i].resize(fMaxNofHitsInEventDraw);
      for (int iH = 0; iH < fMaxNofHitsInEventDraw; iH++){
          fEventHits[i][iH] = new TEllipse(16., 16., 10.);
          fEventHits[i][iH]->SetFillColor(kRed);
          fCirclePicture[i]->AddSpecialObject(fEventHits[i][iH]);
      }
      for (int iH = 0; iH < fMaxNofHitsInEventDraw; iH++){
          fRingHits[i][iH] = new TEllipse(16., 16., 10.);;
          fRingHits[i][iH]->SetFillColor(kYellow+2);
          fCirclePicture[i]->AddSpecialObject(fRingHits[i][iH]);
      }

      AddPicture(fCirclePicture[i], "CirclePicture");
   }

   printf("RICH Histograms created \n");
   fflush ( stdout);
   
   fRingFitter = new CbmRichRingFitterCOP();
   fEllipseFitter = new CbmRichRingFitterEllipseTau();
   fRingFinder = new CbmRichRingFinderHoughImpl();
   
   fEventTime=-1;
   fO2level=-1;
   fH2Olevel=-1;
}

TRICHProc::~TRICHProc()
{
}

void TRICHProc::SetParticleIdentificationConditions()
{
        //21_10_x
        //Double_t electronpntsc1c2[4][2] = { {500, 500}, {450, 1300}, {1200, 1620}, {1300, 500} };
        //30_10_00x
        Double_t electronpntsc1c2[4][2] = { {2150, 720}, {1900, 1800}, {1220, 1730}, {1200, 670} };
        //27_10_017
        //Double_t electronpntsc1c2[4][2] = { {2420, 550}, {2000, 2310}, {1020, 2150}, {940, 540} };
        //23_10_x
        //Double_t electronpntsc1c2[4][2] = { {1375, 500}, {1360, 1540}, {530, 1700}, {520, 500} };
        //29_10_x
        //Double_t electronpntsc1c2[4][2] = { {2180, 645}, {2145, 1920}, {1120, 1910}, {1140, 645} };

        fElectronCondc1c2 = MakePolyCond("RichElectronCondc1c2", 4, electronpntsc1c2);
        fElectronCondc1c2->SetHistogram("Ch1_Ch2");

        Double_t muonpntsc1c2[4][2] = { {500, 500}, {450, 1300}, {1200, 1620}, {1300, 500}  };
        fMuonCondc1c2 = MakePolyCond("RichMuonCondc1c2", 4, muonpntsc1c2);
        fMuonCondc1c2->SetHistogram("Ch1_Ch2");

        Double_t pionpntsc1c2[4][2] = { {500, 500}, {450, 1300}, {1200, 1620}, {1300, 500} };
        fPionCondc1c2 = MakePolyCond("RichPionCondc1c2", 4, pionpntsc1c2);
        fPionCondc1c2->SetHistogram("Ch1_Ch2");
}


bool TRICHProc::ReadMapmtGeometry(const char* fname)
{
  for (int ncell=0; ncell<NUM_MAPMP_CELLS; ncell++) {
    fMAPMTCells[ncell].xbin = 0.;
    fMAPMTCells[ncell].ybin = 0.;
    fMAPMTCells[ncell].xcoord = 0.;
    fMAPMTCells[ncell].ycoord = 0.;
    fMAPMTCells[ncell].globalchannel = 0;
    fMAPMTCells[ncell].mapmt = 0;
    fMAPMTCells[ncell].pixel = 0;
    fMAPMTCells[ncell].roc = 0;
    fMAPMTCells[ncell].feb = 0;
    fMAPMTCells[ncell].febchannel = 0;
  }
  
  std::ifstream f(fname);
  if (!f)  {
    printf("MAPMT geometry file not found %s\n", fname);
    return false;
  }
  
  char sbuf[1024];
  
  // skip first two lines
  f.getline(sbuf, sizeof(sbuf));
  f.getline(sbuf, sizeof(sbuf));
  
  int nread = 0;
  
  while (!f.eof()) {
    f.getline(sbuf, sizeof(sbuf));
    if (strlen(sbuf)==0) break;
    
    int globalchannel(0), mapmt(0), pixel(0), roc(0), feb(0), febchannel(0);
    float xcoord(0.), ycoord(0.), xbin(0.), ybin(0.);
    
    if (sscanf(sbuf, "%d %d %d %d %d %d %f %f %f %f", &globalchannel, &roc, &feb, &febchannel, &mapmt, &pixel, &xcoord, &ycoord, &xbin, &ybin)!=10) {
      printf("MAPMT read line %d format error %s\n", nread, sbuf);
      return false;
    }
    
    int ncell = febchannel;//febchannel + (feb==2 ? 128 : 0);
    
    ncell = roc*256 + febchannel;
    if(feb==2) {
      ncell = ncell + 128;
    }
    
    
    if ((ncell<0) || (ncell>=NUM_MAPMP_CELLS)) {
      printf("MAPMT read wrong cell address %d in line %d\n", ncell, nread);
      return false;
    }
    
    // printf("ncell = %d xbin = %3.1f ybin = %3.1f\n", ncell, xbin, ybin);
    
    fMAPMTCells[ncell].xbin = xbin;
    fMAPMTCells[ncell].ybin = ybin;
    fMAPMTCells[ncell].xcoord = xcoord;
    fMAPMTCells[ncell].ycoord = ycoord;
    fMAPMTCells[ncell].globalchannel = globalchannel;
    fMAPMTCells[ncell].mapmt = mapmt;
    fMAPMTCells[ncell].pixel = pixel;
    fMAPMTCells[ncell].roc = roc;
    fMAPMTCells[ncell].feb = feb;
    fMAPMTCells[ncell].febchannel = febchannel;
    
    nread++;
  }
  
  if (nread!=1024) {
    printf("MAPMT geometry read %d lines, expected %d\n", nread, 1024);
    return false;
  }
  
  printf("Load MAPMT geometry from file %s done\n", fname);
  
  return true;
}

void TRICHProc::FinalizeEvent()
{
  if(!fRocInputEvent->IsValid()) return; // skip incomplete roc events from first step

  // SL: do we need extra processing of pulser events here - probably not
  if (fMbsEvent && fMbsEvent->IsPulser()) return;
  
  ProcessRICH();
  fEventNr++;  
}

void TRICHProc::InitEvent(TGo4EventElement* outevnt)
{
  // first assign input event:
  // since input event object is never discarded within processor lifetime,
  // we just search for subevent by name once to speed up processing  
   if(fRocInputEvent==0) {
      TCBMBeamtimeEvent* btevent=dynamic_cast<TCBMBeamtimeEvent*>(GetInputEvent());
      if(btevent)
      {
         fRocInputEvent=dynamic_cast<TRocEvent*>(btevent->GetSubEvent("ROC"));
         fTriglogEvent=dynamic_cast<TTriglogEvent*>(btevent->GetSubEvent("TRIGLOG"));
         fMbsEvent=dynamic_cast<TMbsCrateEvent*>(btevent->GetSubEvent("MBSCRATE"));
         fEpicsEvent=dynamic_cast<TEpicsEvent*>(btevent->GetSubEvent("EPICS"));
      } 
      else
      {
         fRocInputEvent=dynamic_cast<TRocEvent*>(GetInputEvent());

      }
      if(fRocInputEvent==0) {
         GO4_STOP_ANALYSIS_MESSAGE("**** TRICHProc: Fatal error: input event is not a TRocEvent!!! STOP GO4");
      }
   }

   // then assign output event
   // since output event object is never discarded within processor lifetime,
   // we just search for subevent by name once to speed up processing
   if(fOutputEvent==0)
   {
      TCBMBeamtimeEvent* btevent=dynamic_cast<TCBMBeamtimeEvent*>(outevnt);
      if(btevent)
      {
         fOutputEvent=dynamic_cast<TRICHEvent*>(btevent->GetSubEvent("RICH"));
         fHodo1 = dynamic_cast<TFiberHodEvent*>(btevent->GetSubEvent("Hodo1"));
         fBeamEvent=dynamic_cast<TBeamMonitorEvent*>(btevent->GetSubEvent("BEAM"));
      }
      else
      {

         fOutputEvent= dynamic_cast<TRICHEvent*>(outevnt);
      }
      if(fOutputEvent==0) {
         GO4_STOP_ANALYSIS_MESSAGE("**** TRICHProc: Fatal error: output event is not a TRICHEvent!!! STOP GO4");
      }

   }
}

void TRICHProc::ProcessRICH()
{
   // Get Event time (in secs) from MBS event
   if (fTriglogEvent) {
      fEventTime = fTriglogEvent->fMbsTimeSecs;
   }
   
   // Get O2 and H2O concentration from Epics event, interpolate missing data
   if (fEpicsEvent && fEpicsEvent->IsValid()) {
      fO2level=fEpicsEvent->GetDouble("CBM:RICH:Gas:O2");
      fH2Olevel=fEpicsEvent->GetDouble("CBM:RICH:Gas:H2O");

   }
   if (fEventTime>1318648546 && fEventTime<1318720000) {
      fO2level=350;
      if (fEventTime>1318655338) fO2level=18844646*pow(0.9999834,(fEventTime-1318000000));
   }

  // SaveParameters();
      
   // Get Cherenkov 1/2 signal
   Double_t cherenkov1=fMbsEvent->fData1182[1][0]; //RICH1+2: Electron
   Double_t cherenkov2=fMbsEvent->fData1182[0][1]; //RICH2only: Muon
                                                   //none: Pion
   Double_t leadglass=fMbsEvent->fData1182[1][5];
   
   //perform particle identification
   //Bool_t isElectron = fBeamEvent->fIsElectron;
   //Bool_t isMuon = fBeamEvent->fIsMuon;
   //Bool_t isPion = fBeamEvent->fIsPion;
   bool isElectron = fElectronCondc1c2->Test(cherenkov1,cherenkov2);
   bool isMuon = fMuonCondc1c2->Test(cherenkov1,cherenkov2);
   bool isPion = fPionCondc1c2->Test(cherenkov1,cherenkov2);

   //access information from fiber hodoscope
   Double_t hodx=0;
   Double_t hody=0;
   if (fHodo1 && (fHodo1->NumHits()>0)) {
      hodx = fHodo1->Hit(0).X;
      hody = fHodo1->Hit(0).Y;
   }
   //cout << "hodo x,y " << hodx << " " << hody << endl;
   
   // needed for hodXYscan
   Double_t dx = 64. / fNofHodWindowsX;
   Double_t dy = 64. / fNofHodWindowsY;
   Int_t xNum = (Int_t)( (hodx + 32.) / dx);
   Int_t yNum = (Int_t)( (hody + 32.) / dy);
   
    //printf("\nEventtime in secs: %i \n",fEventTime);
    //printf("O2 content: %f \n",fO2level);
    //printf("H2O content: %f \n",fH2Olevel);
    //printf("Cherenkov 1/2: %f %f \n",cherenkov1,cherenkov2);
    //printf("Hodoscope position: %f / %f \n",hodx,hody);

   Int_t NrHitsperPMT[17];
   for (Int_t iPMT=0; iPMT<=16; iPMT++) NrHitsperPMT[iPMT]=0;
  
   vector<Double_t> timearray;
   //collect all hits from one event in vector
   vector<CbmRichHitLight> hits;
  
   for (UInt_t RICHROC=fPar->RICHFirstROC; RICHROC<=fPar->RICHLastROC; RICHROC++) {
     // cout << "  Processing ROC "<< RICHROC << endl;
      // Loop over all ROCs which deliver RICH data
      
      if (RICHROC>MAX_ROC) continue;
      
      TRocData* RocData=dynamic_cast<TRocData*>(fRocInputEvent->getEventElement(RICHROC));
      if (RocData==0) continue;  //no data from this ROC available;
      
      for (UInt_t MessageNr=0; MessageNr<RocData->fExtMessages.size(); MessageNr++) {
         //cout << "  Processing message " << MessageNr << "from ROC " << RICHROC << endl;
         //cout << "Msg.size() = " << RocData->fExtMessages.size() << endl;
         // Loop over all messages from particular RICH-ROC
         
         TRocMessageExtended msg=RocData->fExtMessages.at(MessageNr);
         if (msg.GetMessageType() != roc::MSG_HIT) continue;  // only interrested in HIT messages
         //cout << "roc::MSG_HIT" << endl;

         // extraxt message information
         UInt_t dataROC=msg.GetRocNumber();
         UInt_t dataFEBNr=msg.GetNxNumber();  // 0 or 2
         UInt_t dataFEBCh=msg.GetNxChNum();   // 0..127
         Double_t dataHitTime=msg.GetTriggerDeltaT();
         Int_t dataHitADCraw=msg.GetNxADC();
         Int_t dataHitADCcor=msg.GetCorrectedNxADC();
                  
         // calculate ncell from ROC and FEB number
         int ncell(0);
         // ncell = dataROC*(RICHROC-fPar->RICHFirstROC) + dataFEBCh;
         ncell=(dataROC- fPar->RICHFirstROC) * 256+dataFEBCh;
         if(dataFEBNr==2) {
                 ncell = ncell + 128;  
         }

         Float_t dist=msg.GetTriggerDeltaT(); // corrected trigger diff is already in input message!
         hHittime->Fill(dist);
         
         // calculate corrected (inverse) ADC value, mapmtPedestal can be changed in TRICHParams.cxx or parameter editor in the Go4 GUI
         int rADC = fPar->mapmtPedestal - msg.GetNxADC();
         
         // fill histogram with all MAPMT ADC values to have reference for condition
         hADCall->Fill(dataHitADCcor);
         
         // get mapping from geometry file
         double xbin = fMAPMTCells[ncell].xbin;
         double ybin = fMAPMTCells[ncell].ybin;
         double xcoord = fMAPMTCells[ncell].xcoord;
         double ycoord = fMAPMTCells[ncell].ycoord;
         Int_t mapmt = fMAPMTCells[ncell].mapmt;
         
         // store mapped info to ROOT tree
         fOutputEvent->fMAPM_X[ncell]=xbin;
         fOutputEvent->fMAPM_Y[ncell]=ybin;
         fOutputEvent->fMAPM_Xcoord[ncell]=xcoord;
         fOutputEvent->fMAPM_Ycoord[ncell]=ycoord;
         fOutputEvent->fMAPM_Integral[ncell]=dataHitADCcor;
         fOutputEvent->fMAPM_Integral_raw[ncell]=rADC;
         fOutputEvent->fMAPM_DeltaT[ncell]=dist;
         
         // printf("GET HIT ncell = %d xbin = %3.1f ybin = %3.1f\n", ncell, xbin, ybin);
         
         // if condition is ok, fill XY histogram.
         // Condition range can be modified in the condition editor in the Go4 GUI
         if (fMAPMT_cond->Test(dataHitADCcor) && fMAPMT_tcond->Test(dist)) {
                 timearray.push_back(dist);

                 // create hit and add it to the vector of hits
                 CbmRichHitLight hit(xcoord, ycoord);
                 hits.push_back(hit);

                 if (xbin < 16 && ybin < 16) hHodXYscan->Fill(xbin + 16. * xNum, ybin + 16. * yNum);
                 
                 if (isElectron) hEntries2d[0]->Fill(xbin, ybin);
                 if (isMuon) hEntries2d[2]->Fill(xbin, ybin);
                 if (isPion) hEntries2d[1]->Fill(xbin, ybin);
                 hEntries2d[3]->Fill(xbin, ybin);

                 Double_t c1=-0.050;
                 Double_t c2=-0.040;
                 Double_t c3=-0.030;
                 Double_t c4=-0.020;
                 Double_t c5=-0.010;
                 hMAPMT_hcor2d_s1->Fill(xbin+c1*hodx,ybin+0*hody);
                 hMAPMT_hcor2d_s2->Fill(xbin+c2*hodx,ybin+0*hody);
                 hMAPMT_hcor2d_s3->Fill(xbin+c3*hodx,ybin+0*hody);
                 hMAPMT_hcor2d_s4->Fill(xbin+c4*hodx,ybin+0*hody);
                 hMAPMT_hcor2d_s5->Fill(xbin+c5*hodx,ybin+0*hody);

                 hEntriesperPMT->Fill(xbin, ybin);
                 hEntries2dIntegral->Fill(xbin, ybin, dataHitADCcor);
                 
                 if (mapmt >=1 && mapmt <= 16) NrHitsperPMT[mapmt]++;
         }
         
         int nbin = hEntries2d[3]->GetBin(hEntries2d[3]->GetXaxis()->FindBin(xbin), hEntries2d[3]->GetYaxis()->FindBin(ybin));
         double entries = hEntries2d[3]->GetBinContent(nbin);
         double integral = hEntries2dIntegral->GetBinContent(nbin);
         hEntries2dMean->SetBinContent(nbin, (entries <= 0.) ? 0. : integral/entries);
      } //MessageNr
   } //RICHROC

   for (Int_t iPMT=1; iPMT<=16; iPMT++) {
      hHitMultperPMT[iPMT]->Fill(NrHitsperPMT[iPMT]);
   }
   
   // time deviation of individual hits from average, not very fast implementation...
   Double_t mean=0;
   for (Int_t ihit=0; ihit<timearray.size(); ihit++){
      // printf("time %i: %f \n",ihit,timearray[ihit]);
      mean=mean+timearray[ihit];
   }
   mean=mean/timearray.size();
   // printf("mean: %f \n",mean);
   for (Int_t ihit=0; ihit<timearray.size(); ihit++) {
      hTimespread->Fill(timearray[ihit]-mean);
   }
   // hit Multiplicity vs O2 level
   if (hits.size() > 5 && hits.size() < 40) hMAPMT_HitMultvsO2->Fill(fO2level, hits.size());

   if (fHodo1 != NULL && fHodo1->NumHits() >=1){
           hNofHitsInHodo->Fill(fHodo1->NumHits());
           hHodoADCvsNofHitInEv->Fill(fHodo1->Hit(0).ADC, hits.size());
   }

   DoAnalysis(hits, isElectron, isPion, isMuon);
}

void TRICHProc::SaveParameters()
{
        ofstream outfile;
        outfile.open ("rich_parameters.txt", ios_base::app);

        double o2 = -1.0;
        double h2o = -1.0;
        double ptb = -1.0;
        double tt1 = -1.0;
        double pt4 = -1.0;
        double refrIndex = -1.0;

        if (fEpicsEvent && fEpicsEvent->IsValid()) {
        o2 = fEpicsEvent->GetDouble("CBM:RICH:Gas:O2");
                h2o = fEpicsEvent->GetDouble("CBM:RICH:Gas:H2O");
        ptb = fEpicsEvent->GetDouble("CBM:RICH:Gas:PTB");
                tt1 = fEpicsEvent->GetDouble("CBM:RICH:Gas:TT-1");
        pt4 = fEpicsEvent->GetDouble("CBM:RICH:Gas:PT-4");
        refrIndex = fEpicsEvent->GetDouble("CBM:RICH:Gas:RefrIndex");
        };

        if (o2 == -1. && h2o == -1. && ptb == -1. && tt1 == -1. && pt4 == -1. && refrIndex == -1.) return;

        TGo4AnalysisStep* firststep = TGo4Analysis::Instance()->GetAnalysisStepNum(0);
        if (dynamic_cast<TGo4MbsFileParameter*>(firststep->GetEventSource())) {
                outfile << firststep->GetEventSourceName() << " " <<
                                fEventTime << " " <<
                                o2 << " " <<
                                h2o << " " <<
                                ptb << " " <<
                                tt1 << " " <<
                                pt4 << " " <<
                                refrIndex << endl;
        }
        outfile.close();
        time_t rawtime = fEventTime;
        //time ( &rawtime );
        printf ( "The current local time is: %s", ctime (&rawtime) );
}

void TRICHProc::DoAnalysis(
                const vector<CbmRichHitLight>& hits,
                bool isElectron,
                bool isPion,
                bool isMuon)
{
        //cout << "hits.size() = " << hits.size() << endl;
        if (hits.size() < 1) return;

    //[0] is electron, [1] is pion, [2] is muon, [3] is all
        Int_t hIndex = -1;
        if (isElectron) hIndex = 0;
        else if (isPion) hIndex = 1;
        else if (isMuon) hIndex = 2;
        if (hIndex < 0) hIndex = 2;

        hNofHitsEvAll[hIndex]->Fill(hits.size());
        hNofHitsEvAll[3]->Fill(hits.size());

        fRingFinder->DoFind(hits);
        vector<CbmRichRingLight*> rings = fRingFinder->GetFoundRings();
        hNofFoundRings->Fill(rings.size());
        if (rings.size() <= 0) return;
        CbmRichRingLight* ring = rings[0];
        //if (rings.size() == 2) cout << "2 rings were found" << endl;

        Int_t nofHits = ring->GetNofHits();
        for (Int_t iH = 0; iH < hits.size(); iH++){
                Bool_t hitFound = false;
                for (Int_t i = 0; i < nofHits; i++) {
                   if (ring != NULL && ring->GetHit(i).fX == hits[iH].fX && ring->GetHit(i).fY == hits[iH].fY){
                           hitFound = true;
                   }
                }
                if (hitFound == false) hNotAssignedHitsXY->Fill(hits[iH].fX, hits[iH].fY);
        }

        hNofHitsRingFinder[hIndex]->Fill(ring->GetNofHits());
        hNofHitsRingFinder[3]->Fill(ring->GetNofHits());

        hNofHitsEvRingFinder[hIndex]->Fill(hits.size());
        hNofHitsEvRingFinder[3]->Fill(hits.size());
        hNofHitsEvRingFinderEff[hIndex]->Divide(hNofHitsRingFinder[hIndex], hNofHitsEvAll[hIndex]);
        hNofHitsEvRingFinderEff[hIndex]->Divide(hNofHitsRingFinder[3], hNofHitsEvAll[3]);
        //cout << "Ring finder eff:" << 100.* (double)hNofHitsRingFinder[3]->GetEntries()/ hNofHitsAll[3]->GetEntries() << endl;

        //fit COP
        for (int iR = 0; iR < rings.size(); iR++){
                fRingFitter->DoFit(rings[iR]);
        }
        FillCircleFitHistograms(hIndex, ring, hits.size());// fill histograms for a specific particle
        FillCircleFitHistograms(3, ring, hits.size());// fill histograms for all particle. [3] is all
        DrawEventHits(hits, fNrSED);
        DrawRingHits(ring, fNrSED);
        DrawCircles(rings, fNrSED);

        for (int iR = 0; iR < rings.size(); iR++){
                fEllipseFitter->DoFit(rings[iR]);
        }
        if (hIndex>=0) FillEllipseFitHistograms(hIndex, ring, hits.size());// fill histograms for a specific particle
        FillEllipseFitHistograms(3, ring, hits.size());// fill histograms for a specific particle
        DrawEllipses(rings, fNrSED);

        // rotate single event printing to next histogram ONLY if last one had a non-empty event
        if ( rings.size() == 2 ){// && ring->GetBaxis()/ring->GetAaxis() < 0.7) {
                fNrSED++;
                if (fNrSED >= hSED.size()) fNrSED = 0;
        }
        //if (ring != NULL) delete ring;
}

void TRICHProc::FillCircleFitHistograms(
                Int_t hIndex,
                CbmRichRingLight* ring,
                int nofEventHits)
{
   if (hIndex != -1){
      if (ring->GetRadius() > 0. && ring->GetCenterX() > 0. && ring->GetCenterY() > 0.) {
         hCircleFitRadius[hIndex]->Fill(ring->GetRadius());
         hCircleFitCenter[hIndex]->Fill(ring->GetCenterX(), ring->GetCenterY());
         hCircleFitChi2[hIndex]->Fill(ring->GetChi2() / ring->GetNofHits());
         for (int iH = 0; iH < ring->GetNofHits(); iH++){
                 double xc = ring->GetCenterX();
                 double yc = ring->GetCenterY();
                 CbmRichHitLight h = ring->GetHit(iH);
                 double dr = sqrt((xc - h.fX)*(xc - h.fX) + (yc - h.fY)*(yc - h.fY)) - ring->GetRadius();
                 hCircleFitDR[hIndex]->Fill(dr);
         }
      }
      // good fitted ring
      if (ring->GetRadius() > 20 && ring->GetRadius() < 70) hNofHitsEvCircleFit[hIndex]->Fill(nofEventHits);
      hNofHitsEvCircleFitEff[hIndex]->Divide(hNofHitsEvCircleFit[hIndex], hNofHitsEvAll[hIndex]);
   }
}

void TRICHProc::FillEllipseFitHistograms(
                Int_t hIndex,
                CbmRichRingLight* ring,
                int nofEventHits)
{
   if (ring->GetAaxis() > 0. && ring->GetBaxis() > 0. && ring->GetPhi() != -1. && ring->GetCenterX() != -1. && ring->GetCenterY()!=-1.) {
      hEllipseFitAaxis[hIndex]->Fill(ring->GetAaxis());
      hEllipseFitBaxis[hIndex]->Fill(ring->GetBaxis()); 
      hEllipseFitBoverA[hIndex]->Fill(ring->GetBaxis()/ring->GetAaxis());  
      hEllipseFitPhi[hIndex]->Fill(ring->GetPhi());    
      hEllipseFitCenter[hIndex]->Fill(ring->GetCenterX(),ring->GetCenterY());
      hEllipseFitChi2[hIndex]->Fill(ring->GetChi2() / ring->GetNofHits());
   }
   //good fitted ring with ellipse fitter
   if (ring->GetAaxis() > 20 && ring->GetAaxis() < 80 && ring->GetBaxis() > 20
                   && ring->GetBaxis() < 80) hNofHitsEvEllipseFit[hIndex]->Fill(nofEventHits);
   hNofHitsEvEllipseFitEff[hIndex]->Divide(hNofHitsEvEllipseFit[hIndex], hNofHitsEvAll[hIndex]);
}

void TRICHProc::DrawRingHits(
                CbmRichRingLight* ring,
                Int_t sedNum)
{
        if (ring->GetNofHits() > fMaxNofHitsInEventDraw) return;

   Int_t nofHits = ring->GetNofHits();
   for (Int_t i = 0; i < nofHits; i++) {
          fRingHits[sedNum][i]->SetX1(ring->GetHit(i).fX);
          fRingHits[sedNum][i]->SetY1(ring->GetHit(i).fY);
          fRingHits[sedNum][i]->SetR1(3);
          fRingHits[sedNum][i]->SetR2(3);
   }
   for (Int_t i = nofHits; i < fMaxNofHitsInEventDraw; i++){
           fRingHits[sedNum][i]->SetX1(-999999);
   }
}

void TRICHProc::DrawEventHits(
                const vector<CbmRichHitLight>& hits,
                Int_t sedNum)
{
   if (hits.size() > fMaxNofHitsInEventDraw) return;
   Int_t nofHits = hits.size();
   for (Int_t i = 0; i < nofHits; i++) {
           fEventHits[sedNum][i]->SetX1(hits[i].fX);
           fEventHits[sedNum][i]->SetY1(hits[i].fY);
           fEventHits[sedNum][i]->SetR1(4);
           fEventHits[sedNum][i]->SetR2(4);
   }
   for (Int_t i = nofHits; i < fMaxNofHitsInEventDraw; i++){
           fEventHits[sedNum][i]->SetX1(-999999);
   }
}

void TRICHProc::DrawCircles(
                const vector<CbmRichRingLight*>& rings,
                Int_t sedNum)
{
        if (rings.size() <= 0) return;
        for (int iR = 0; iR < rings.size(); iR++){
           fCircle[sedNum][iR]->SetX1(rings[iR]->GetCenterX());
           fCircle[sedNum][iR]->SetY1(rings[iR]->GetCenterY());
           fCircle[sedNum][iR]->SetR1(rings[iR]->GetRadius());
           fCircle[sedNum][iR]->SetR2(rings[iR]->GetRadius());
   }
   stringstream ss;
   ss.precision(3);
   ss << "(x,y,r,n):(" << rings[0]->GetCenterX() << ", " << rings[0]->GetCenterY() << ","
                   << rings[0]->GetRadius()<< ", " << rings[0]->GetNofHits() << ")";
   fCircleText[sedNum]->SetText(5, 200, ss.str().c_str());
}


void TRICHProc::DrawEllipses(
                const vector<CbmRichRingLight*>& rings,
                Int_t sedNum)
{
        for (int iR = 0; iR < rings.size(); iR++){
           fEllipse[sedNum][iR]->SetX1(rings[iR]->GetCenterX());
           fEllipse[sedNum][iR]->SetY1(rings[iR]->GetCenterY());
           fEllipse[sedNum][iR]->SetR1(rings[iR]->GetAaxis());
           fEllipse[sedNum][iR]->SetR2(rings[iR]->GetBaxis());
           fEllipse[sedNum][iR]->SetTheta(rings[iR]->GetPhi() * 180. / 3.1415);
        }
   stringstream ss;
   ss.precision(3);
   ss << "(x,y,a,b,p,n):(" << rings[0]->GetCenterX() << ", " << rings[0]->GetCenterY() << ","
      << rings[0]->GetAaxis()<< ", "<< rings[0]->GetBaxis()<< ", " << rings[0]->GetPhi()<< ", "
      << rings[0]->GetNofHits() << ")";
   fEllipseText[sedNum]->SetText(5, 185, ss.str().c_str());
}