HarpoAnalyseRates.cxx

Go to the documentation of this file.

// 
// File HarpoAnalyseRates.cxx
//
#include "HarpoAnalyseRates.h"
#include "HarpoConfig.h"
#include "HarpoDetSet.h"
#include "HarpoDebug.h"
#include "HarpoDccEvent.h"
#include "Pmm2Event.h"
#include "Pmm2Header.h"
#include "HarpoEvent.h"
#include "HarpoRecoEvent.h"

#include "TFile.h"
#include "TStyle.h"
#include "TCanvas.h"
#include "TLatex.h"
#include "TGraph.h"
#include "TF1.h"
#include "TMath.h"
#include "TRandom.h"
#include "TSystem.h"

#include <cstdlib>
#include <cstring>
#include <cassert>
#include <fstream>
#include <iostream>

ClassImp(HarpoAnalyseRates);

void   HarpoAnalyseRates::print()
 {
   unsigned int nd; // number of detectors
   HarpoEventHeader *hdr = fEvt->GetHeader();

   assert(hdr != NULL);
   hdr->print();
   
   for (nd = 0; nd < gkNDetectors; nd++) {
     //     if (fEvt->isdataExist(nd)) {
       HarpoDccMap *plane = fEvt->GetDccMap(nd);
       if (plane != NULL )plane->print();
   }
 }

void   HarpoAnalyseRates::process()
{
  //  Bool_t drawEvent = kFALSE;
  nEvents++;
  
  if(!gHDetSet->isExist(PMM2))
    return;

  if(nEvents%10000 == 0)
    Info("process","Event %ld", nEvents);

  Pmm2Event *anaEvt = static_cast<Pmm2Event *>(fEvt->GetDetEvent(PMM2));
 
  Double_t timestamp = fEvt->GetTimeStamp();
  if(!timestamp) timestamp = 10*anaEvt->GetTimeStamp();
  if(timestamp==0)
    return;
  //  Info("process","TS %g", timestamp);

  if(fInit == kFALSE){
    fInit = kTRUE;

    fTimeStart = timestamp;
    
    Pmm2Header* pmm2Hdr = (Pmm2Header*)fRunHeader->GetHeader(PMM2);
    fTriggerTimeTot = 0;
    for(Int_t i = 0; i<kNtriggerTypes; i++){
      fTriggerTime[i] = 0;
      fTriggerDownscaling[i] = 1;
      if(pmm2Hdr == NULL) continue;
      Pmm2Status* pmm2Stat = pmm2Hdr->getStatus();
      if(pmm2Stat == NULL) continue;
      unsigned int triggerinfo = pmm2Stat->getTrigger(i); // 32bits: 
      fTriggerTime[i] = triggerinfo/65536;
      fTriggerDownscaling[i] = triggerinfo%65536;
    }
    for(Int_t i = 0; i<kNtriggerTypes; i++)
      fTriggerTimeTot += fTriggerTime[i];
    Info("Init","TriggerTimeTot = %d", fTriggerTimeTot);
    for(Int_t i = 0; i<kNtriggerTypes; i++){
      Double_t trigCorr  = 1.;
      if(fTriggerTime[i])
        trigCorr /= fTriggerTime[i]*1./fTriggerTimeTot;
      if(fTriggerDownscaling[i])
        trigCorr *= fTriggerDownscaling[i];
      hTriggerCorrection->SetBinContent(i+3,trigCorr);
      Info("Init","Trigger param %d %d => %g",fTriggerTime[i],fTriggerDownscaling[i],trigCorr);
    }
    hTriggerCorrection->SetBinContent(1,1);
    hTriggerCorrection->SetBinContent(2,1);
    hTriggerCorrection->SetBinContent(3,1);
    hCycle = new TH2F("hCycle",";t;trigger",fTriggerTimeTot,0,fTriggerTimeTot*100000,20,-2,18);
    hCycleNew = new TH2F("hCycleNew",";t;trigger",fTriggerTimeTot,0,fTriggerTimeTot*100000,20,-2,18);
    //hTlowhigh = new TH2F("hTlowhigh",";t_{low};t_{high}",fTriggerTimeTot,0,fTriggerTimeTot*100000,fTriggerTimeTot,0,fTriggerTimeTot*100000);
    for(Int_t i = 0; i<kNtriggerTypes; i++){
      fTlow[i] =  fTriggerTimeTot*100000;
      fThigh[i] = 0;
    }
  }
  
  //timestamp -= fTimeStart;


  
  Int_t triggerType = -1;

  // Pmm2Event *anaEvt = static_cast<Pmm2Event *>(fEvt->GetDetEvent(PMM2));
  triggerType = anaEvt->GetTriggerType(); // Trigger line (with mesh, without, traversing tracks, ...)
  
  int imes;
  int esize = anaEvt->GetHeader()->eventSize;   // number of PMs hit
  Pmm2MesVect * pm = anaEvt->GetMesurements();  // data for each included PM
  for(imes=0;imes<esize;imes++) {
    //    int tsPm = pm->at(imes).getTimeStamp(); // charge
    Double_t tsPm = pm->at(imes).getTimeStamp() + 16777216*Double_t(pm->at(imes).getPPS());
    hTimestampPmm2->Fill(Double_t(tsPm)*10. - Double_t(timestamp+fTimeStart));
    //Info("process","%d - %d/10 = %g",Int_t(tsPm),Int_t(timestamp),Double_t(tsPm) - Double_t(timestamp)/10.);
  } 

  hCycle->Fill(fmod(timestamp,fTriggerTimeTot*100000.),triggerType);

  // if(triggerType == -1)
  //   Info("process","Unknown trigger %d %d",triggerType,fTriggerOld);
  //  if(triggerType != fTriggerOld && fTriggerOld != -1){
  if(triggerType == -1){
    Int_t tLow = fTlow[fTriggerOld];
    // Int_t iMax = fTriggerOld - 8;
    // if(iMax<0) iMax += kNtriggerTypes;
    // for(Int_t i = 0; i<=iMax; i++){
    //  tLow += fTriggerTime[i+8]*100000;
    //  if(tLow>fTriggerTimeTot*100000) tLow -= fTriggerTimeTot*100000;
    // }
    //Info("process","Unknown trigger => correction %.3g %.3g %.3g",fmod(timestamp,fTriggerTimeTot*100000.)*1.,tLow*1.,fmod(timestamp,fTriggerTimeTot*100000.)*1./tLow);
    if(fmod(timestamp,fTriggerTimeTot*100000.) < tLow){
      triggerType = fTriggerOld;
      hChangeTrigger->Fill(fTriggerOld);
    }
    for(Int_t i = 0; i<kNtriggerTypes; i++){
      if(fmod(timestamp,fTriggerTimeTot*100000.) > fTlow[i] && fmod(timestamp,fTriggerTimeTot*100000.) < fThigh[i]){
        triggerType = i;
        hChangeTrigger->Fill(i);
      }
    }
  }
  // }else{

  if( (triggerType-fTriggerOld)==1 || (triggerType-fTriggerOld)==kNtriggerTypes-1){
    Int_t tLow = fmod(fTimestampOld,fTriggerTimeTot*100000.);
    Int_t tHigh = fmod(timestamp,fTriggerTimeTot*100000.);
    //Info("process","**** %d %d",tLow,tHigh);
    // Int_t iMax = triggerType - 8;
    // if(iMax<0) iMax += kNtriggerTypes;
    // for(Int_t i = 0; i<iMax; i++){
    //   tLow -= fTriggerTime[i+8]*100000;
    //   tHigh -= fTriggerTime[i+8]*100000;
    //   if(tLow<0) tLow += fTriggerTimeTot*100000;
    //   if(tHigh<0) tHigh += fTriggerTimeTot*100000;
    //   //Info("process","     %d %d",tLow,tHigh);
    // }
    if(tLow<fTlow[triggerType]) fTlow[triggerType] = tLow;
    if(tHigh>fThigh[fTriggerOld]) fThigh[fTriggerOld] = tHigh;
    //      hTlowhigh->Fill(tLow,tHigh);
    // gTlow->SetPoint(gTlow->GetN(),timestamp,fTlow[fTriggerOld]);
    // gThigh->SetPoint(gThigh->GetN(),timestamp,fThigh[fTriggerOld]);
    gTlow->SetPoint(gTlow->GetN(),fTlow[fTriggerOld],fTriggerOld + timestamp*1e-12);
    gThigh->SetPoint(gThigh->GetN(),fThigh[fTriggerOld],fTriggerOld + timestamp*1e-12);
  }
   //Info("process","Trigger change %d %d, %d %d",triggerType,fTriggerOld,fTlow,fThigh);
  //  }

  hCycleNew->Fill(fmod(timestamp,fTriggerTimeTot*100000.),triggerType);


  if(triggerType == fTriggerOld || triggerType == -1)
    hDeltaTlog->Fill(triggerType,Double_t(timestamp - fTimestampOld)/1e5 + 0.05*gRandom->Rndm());
  if(triggerType == -1){
    triggerType = fTriggerOld;

    // for(Int_t i = 0; i<kNtriggerTypes; i++)
    //   fTriggerTimeTot += fTriggerTime[i];

  }

  if(triggerType == fTriggerOld){
    hDeltaT->Fill(triggerType,Double_t(timestamp - fTimestampOld)/1e5);
  }
  else{
    hDeltaTchange->Fill(triggerType,fTriggerOld,Double_t(timestamp - fTimestampOld)/1e5);
    hTriggerChange->Fill(triggerType,fTriggerOld);
  }




  hRate->Fill(triggerType);

  if(triggerType == 0 && fTriggerOld != 0)
    fOffset = timestamp;

  //  hCycle->Fill(timestamp-fOffset,triggerType);

  // gTimestamp0->SetPoint(gTimestamp0->GetN(),nEvents,fEvt->GetTimeStamp(XDCC));
  // //  gTimestamp1->SetPoint(gTimestamp1->GetN(),nEvents,fEvt->GetTimeStamp(YDCC));
  // gTimestamp1->SetPoint(gTimestamp1->GetN(),nEvents,timestamp);
  // gTimestamp2->SetPoint(gTimestamp2->GetN(),nEvents,fEvt->GetTimeStamp(PMM2));

  fTimestampOld = timestamp;
  if(triggerType >=0)
    fTriggerOld = triggerType;

}

void   HarpoAnalyseRates::Init()
{

  fInit = kFALSE;

  fTimestampOld = 0;
  fTriggerOld = -4;
  fOffset = 0;
  fTimeStart = 0;

  const Int_t nbinsDt = 1000;
  Double_t xminDt = 1;
  Double_t xmaxDt = 1e4;
  Double_t logxminDt = TMath::Log(xminDt);
  Double_t logxmaxDt = TMath::Log(xmaxDt);
  Double_t binwidthDt = (logxmaxDt-logxminDt)/nbinsDt;
  Double_t xbinsDt[nbinsDt+1];
  xbinsDt[0] = xminDt;
  for (Int_t i=1;i<=nbinsDt;i++)
    xbinsDt[i] = TMath::Exp(logxminDt+i*binwidthDt);


  // Initialise histograms here
  hTimestampPmm2 = new TH1F("hTimestampPmm2","t_{PMM2} - t",50,-1310,-1260);
  hDeltaT = new TH2F("hDeltaT",";trigger;#Deltat [ms]",20,-2,18,1000,0,500);
  hDeltaTlog = new TH2F("hDeltaTlog",";trigger;#Deltat [ms]",20,-2,18,nbinsDt,xbinsDt);
  hRate = new TH1F("hRate",";trigger;Rate",20,-2,18);
  hTriggerCorrection = new TH1F("hTriggerCorrection",";trigger;Rate",20,-2,18);
  for(Int_t i = 1; i<= 20; i++)
    hTriggerCorrection->SetBinContent(i,1);
  hDeltaTchange = new TH3F("hDeltaTchange","trigger_{n};trigger_{n-1};#Deltat [ms]",20,-2,18,20,-2,18,1000,0,100);
  hTriggerChange = new TH2F("hTriggerChange",";trigger_{n};trigger_{n-1};",20,-2,18,20,-2,18);
  hChangeTrigger = new TH1F("hChangeTrigger",";changed trigger;",20,-2,18);
  


  gTimestamp0 = new TGraph();
  gTimestamp1 = new TGraph();
  gTimestamp2 = new TGraph();

  gTlow = new TGraph();
  gThigh = new TGraph();




 }

void   HarpoAnalyseRates::Save(char * /* mode */)
 {



   hTriggerCorrection->Multiply(hRate);
   Double_t scale = 1.;
   if(fTimestampOld)
     scale = 1e8/fTimestampOld;
   Info("Save","scale = 1e8/%g = %g",fTimestampOld,scale);
   hTriggerCorrection->Scale(scale);
   Info("Save","fTimestampOld = %g",fTimestampOld);


   /* TFile* hf = */ OpenHistFile("rates");
   // Save histograms here
   hDeltaT->Write();
   hDeltaTlog->Write();
   hDeltaTchange->Write();
   hTriggerChange->Write();
   hTimestampPmm2->Write();
   if(fInit){
     hCycle->Write();
     hCycleNew->Write();
     //     hTlowhigh->Write();
   }
   hChangeTrigger->Write();
   hRate->Write();
   hTriggerCorrection->Write();
   gTimestamp0->Write("gTimestamp0");
   gTimestamp1->Write("gTimestamp1");
   gTimestamp2->Write("gTimestamp2");
   gTlow->Write("gTlow");
   gThigh->Write("gThigh");
 
 }