mcr/calib/tofino/CountNumberOfEvents.C

#include "TH2.h"
#include "TCanvas.h"
#include "TFile.h"
#include "TApplication.h"
#include "TTree.h"
#include "TF1.h"
#include "TNamed.h"
#include "TGraph.h"
#include "TProfile.h"

#include "TFndHdt.h"
//#include "FROOT.h"

using namespace FROOT;

TFile *fHdtFile = 0;

TFndHdt *fCurHdtEv = 0;
TTree *fHdtTree = 0;

TProfile *Nbhabha[K_N_TOFINO_SLABS], *Nkaon[K_N_TOFINO_SLABS];
Int_t N_bhabha[K_N_TOFINO_SLABS], N_kaon[K_N_TOFINO_SLABS];

void DeclareHistos(Int_t run_number_start, Int_t run_number_stop){

  TString hname;
  TString htitle;

  Int_t N=run_number_stop-run_number_start+1;

  for (Int_t sl=0; sl<12; sl++){
    hname.Form("Nbhabha_%d",sl+1);
    htitle.Form("Number of bhabha events TOFINO slab %d", sl+1);
    Nbhabha[sl] = new TProfile(hname, htitle, N, run_number_start-0.5, run_number_stop+0.5);
    
    hname.Form("Nkaon_%d",sl+1);
    htitle.Form("Number of K+K- events TOFINO slab %d", sl+1);
    Nkaon[sl] = new TProfile(hname, htitle, N, run_number_start-0.5, run_number_stop+0.5);
  }
}

void CloseHdtFile(){
  
  if (fHdtTree){
    delete fHdtTree;
    // fHdtTree->Delete();
    fHdtTree = 0;
  }
  
  if(fHdtFile){
    if(fHdtFile->IsOpen()) fHdtFile->Close();
    delete fHdtFile;
    fHdtFile = 0;
  }
  
}

Int_t OpenHdtFile(const TString &runtype,const Int_t &run_num){
  // return value:
  //              0: ok
  //             -1: file not found

  CloseHdtFile();

  TString runname = BuildRunName(runtype,run_num);
  TString fnam;
  fnam.Form("%s/%s.hdt.root",ExpandPathName("$FND_HDT").Data(),runname.Data());
  
  if(gSystem->AccessPathName(fnam)){
    gROOT->Warning("OpenHdtFile","file \"%s\" not found",fnam.Data());
    return -1;
  }

  gROOT->Info("OpenHdtFile","opening file \"%s\"",fnam.Data());
  fHdtFile = new TFile(fnam,"OPEN");

  return 0;
}


Int_t LoadHdtTree(){
  // return value:
  //              0: ok
  //             -1: Problems while getting HDT tree from file
  
  if (fHdtTree){
    delete fHdtTree;
    // fHdtTree->Delete();
    fHdtTree = 0;
  }
  
  fHdtTree = (TTree *) ( fHdtFile->Get("FIN_HDT_TREE") );
  if(!fHdtTree){
    gROOT->Warning("LoadHdtTree","Problems while getting HDT tree from file");
    return -1;
  }
  //  fHdtTree->Print();

  fHdtTree->SetBranchAddress("fndhdt",&fCurHdtEv);

  return 0;
}

Short_t GetTrigType(){
  TFndTrig *cur_trig = (TFndTrig *) (fCurHdtEv->GetGts());
  Short_t hyp_bit=cur_trig->GetPU_cont_bit(0, 0);
  Short_t bha_bit=cur_trig->GetPU_cont_bit(0, 1);
  return hyp_bit*10+bha_bit;
}

void GetMultiplicity(Int_t *mult){
  // mult[0] ... low threshold
  // mult[1] ... high threshold
  mult[0]=0; mult[1]=0;

  TClonesArray *TofHits = fCurHdtEv->GetTofHits();
  Int_t lay, slab;
  TFndHTof *cur_hit = 0;

  for(Int_t ic=0;ic< TofHits->GetEntries(); ic++) {
    cur_hit =  (TFndHTof *) ( TofHits->At(ic) );
    cur_hit->GetLayer_Slab(lay, slab);
    if(lay==E_FIN_OUTER_LAYER) continue;

    Bool_t both_adc = kFALSE;
    Bool_t both_tdcl = kFALSE;
    Bool_t both_tdch = kFALSE;
    if(cur_hit->GetBothAdc() == TFndHit::E_HIT_BOTH_CH ) both_adc = kTRUE;
    if(cur_hit->GetBothLowTdc() == TFndHit::E_HIT_BOTH_CH ) both_tdcl = kTRUE;
    if(cur_hit->GetBothHighTdc() == TFndHit::E_HIT_BOTH_CH ) both_tdch = kTRUE;

    Double_t adce=cur_hit->GetAdcE();
    Double_t adcp=cur_hit->GetAdcP();

    if(both_tdch&&adce>0&&adcp>0){ // kaons
      mult[1]++;
    }   
    if(both_tdcl&&adce>0&&adcp>0){ // bhabha+kaons
      mult[0]++;
    }
  }
}
  
void FillHistos(const Int_t &run_num){
  
  TClonesArray *TofHits = fCurHdtEv->GetTofHits();
  
  Int_t lay, slab;
  TFndHTof *cur_hit = 0;
  Int_t mult[2];
  GetMultiplicity(mult);
  for(Int_t ic=0;ic< TofHits->GetEntries(); ic++) {
    cur_hit =  (TFndHTof *) ( TofHits->At(ic) );
    //    cur_hit->PrintHit();
    cur_hit->GetLayer_Slab(lay, slab);
    //     TThread::Printf("%d %d %d",lay, cham, wire);
    //     TThread::Printf("");
    if(lay==E_FIN_OUTER_LAYER) continue;

    Bool_t both_adc = kFALSE;
    Bool_t both_tdcl = kFALSE;
    Bool_t both_tdch = kFALSE;
    if(cur_hit->GetBothAdc() == TFndHit::E_HIT_BOTH_CH ) both_adc = kTRUE;
    if(cur_hit->GetBothLowTdc() == TFndHit::E_HIT_BOTH_CH ) both_tdcl = kTRUE;
    if(cur_hit->GetBothHighTdc() == TFndHit::E_HIT_BOTH_CH ) both_tdch = kTRUE;

    //     cout << "a: " << cur_hit->GetBothAdc() << " ; l: " << cur_hit->GetBothLowTdc() << " ; h: " << cur_hit->GetBothHighTdc() << endl;
    //     if(both_tdcl) continue; // pedestal
    //     if(!both_tdcl || both_tdch) continue; // bhabha
    //     if(!both_tdcl) continue; // bhabha + kaons
    //    if(!both_tdch) continue; // kaons

    Double_t adce=cur_hit->GetAdcE();
    Double_t adcp=cur_hit->GetAdcP();

    if(both_tdcl && !both_tdch){ // bhabha 
      if (GetTrigType()==1&&mult[0]==2) N_bhabha[slab-1]++;
    }   

    if(both_tdch){ // kaons
      if (GetTrigType()==10&&mult[1]==2) N_kaon[slab-1]++;
    }
  }
}

void EventLoop(const Int_t &run_num, const Int_t &nevs){
  
  Int_t entries = fHdtTree->GetEntries();
  
  Int_t cur_ev = 0;
  while(1){
    if(cur_ev >= entries){      
      gROOT->Info("EventLoop","run completed");
      break;
    }
    if(nevs != -1 && cur_ev >= nevs){      
      gROOT->Info("EventLoop","requested number of events reached");
      break;
    }
    //
    fHdtTree->GetEntry(cur_ev);
    //
    FillHistos(run_num);
    
    if(cur_ev%3000 == 0 ) cout << "cur_ev: " << cur_ev << endl;
    //     fCurHdtEv->PrintTofHits();
    //
    cur_ev++;
  }
}

void SaveHistos(const TString &fname){
  TFile *fOutFile = new TFile(fname,"RECREATE");

  for(Int_t sl=0;sl<K_N_TOFINO_SLABS;sl++ ){
    Nbhabha[sl]->Write();
    Nkaon[sl]->Write();
   }
  
//   TofCalibration->Write();
//   cout << "================> " << TofCalibration->GetPedE(1) << endl;
//   cout << "================> " << TofCalibration->GetPedE(4) << endl;

  fOutFile->Close();
  delete fOutFile;
  fOutFile = 0;

}


Int_t GetRunInfo(TString runtype,Int_t runnum){

  TString run_path;
  run_path.Form("%s/%s.hdt.root",FROOT::ExpandPathName("$FND_HDT").Data(),FROOT::BuildRunName(runtype,runnum).Data());

  TFile *hdt_file = new TFile(run_path);
  if (!hdt_file) return -1;

  //  hdt_file->ls();
  TFndGenInfo *geninfo = (TFndGenInfo *)hdt_file->Get("TFndGenInfo");
  
  //  geninfo->Print();
  
  TDatime d = geninfo->GetDateTime();
  
  return (Int_t) d.Convert();

}

void CountNumberOfEvents(const TString &runtype,const Int_t &run_num_start,const Int_t &run_num_stop,const TString &root_filename){
  
  DeclareHistos(run_num_start, run_num_stop);
  // --- loop in runs
  for(Int_t cur_run=run_num_start; cur_run <=run_num_stop ; cur_run++){
    cout << "Filling from run \"" << runtype.Data() << "\" " << cur_run << endl;
    for (Int_t sl=0; sl<K_N_TOFINO_SLABS; sl++){
    N_bhabha[sl]=0; N_kaon[sl]=0;
    }

    if(OpenHdtFile(runtype,cur_run) !=0) continue; 
    if(LoadHdtTree() != 0) continue;

    EventLoop(cur_run, -1);      
    CloseHdtFile();
    for (Int_t sl=0; sl<K_N_TOFINO_SLABS; sl++){
      Nbhabha[sl]->Fill(cur_run, N_bhabha[sl]);
      Nkaon[sl]->Fill(cur_run, N_kaon[sl]);
    }
  }
  SaveHistos(root_filename);
}

---