PREAN/TFndProcessRec.cxx

// @(#)fROOT/PREAN:$Name:  $:$Id: TFndProcessRec.cxx,v 1.18 2007/09/24 07:32:41 Diego_Faso Exp $
// Final Revision Author: Diego Faso <mailto:faso@to.infn.it>, 2006/10/30
// First Revision Author: B. Dalena, D. Di Santo, D. Faso, G. Simonetti
// Original Author: Valerio Filippini, Alberto Panzarasa (2000)

//                                                                                        //
//    TFndProcessRec:                                                                     // 
//      root interface to fortran FINUDA reconstruction program (fidarc)                  //
//                                                                                        //
//    06.12.2003: Class updated and modified in turin for fROOT by:                       //
//                Barbara Dalena, Daniela Di Santo, Diego Faso, Giuseppe Simonetti        //
//                                                                                        //
// Every change in the reconstruction proceduremust be applied here if some of the        //
// following routines are changed:                                                        //
//  - main.f -----> check "InitRun()"                                                     //
//  - uginit.f ---> just check fortran (if something new has been added)                  //
//  - gutrev.f ---> check "FillJFGES()" up to "IF (HYPE) THEN" [in gutrev.f]              //
//             ---> check "Reconstruct()" from "IF (HYPE) THEN" [in gutrev.f]             //
//                  for kaons: remember to update both cases: "HYPE" and "HYPE && BHABHA" //
//                                                                                        //
// zebra access example methods                                                           //
// // --- you can test the class like this:                                               //
//                                                                                        //
//    TFndProcessRec *proc =new TFndProcessRec("FINU",536,kTRUE);                         //
//    proc->SetMaxEv(100);                                                                //
//    proc->AutoProcessRun();                                                             //
//                                                                                        //
// Automatic file 'Autoprocess_histos.root' is created                                    //
// BHABHA bank not decoded in this example!!!                                             //

#include <dlfcn.h>
#include <Riostream.h>
#include <TROOT.h>
#include <TSystem.h>
#include <TMath.h>

#include "TFndProcessRec.h"

// --- Used for examples
#include "TCanvas.h"
#include "TFile.h"
// ---
const Int_t Ng=1500000; // good
const Int_t Nh=500000; // good
// from main.f
//  (NG = 5000024, NH = 1000000) // before fida-version 603
//  (NG = 5000024, NH = 5000024) // from fida-version 603

struct FidaIntPawc   { int H[Nh]; } pawc_;
struct FidaIntGcBank { int Q[Ng]; } gcbank_;

ClassImp(TFndProcessRec)

//_______________________________________________________________
TFndProcessRec::TFndProcessRec()
:TNamed("dummyinterface","dummy interface"),
 fProcDirBaseName(),fCurProcNum(),fTyp(),fNum(),fRawFileName(),fMsqlDbHostStr(),
 fZebini(),fGcflag(),fRuninf(),flib_handle(),
 fJINFO(),fJFGES(),fJFDST(),fLQ(),fIQ(),fQ(),fdebug(),fRawFileExists()
{
  //
  CommInit();
}

//_______________________________________________________________
TFndProcessRec::TFndProcessRec(Int_t num, Bool_t debug)
:TNamed("froot<->fidarc_interface","FINUDA froot<->mc_interface"),
 fProcDirBaseName(),fCurProcNum(),fTyp(),fNum(),fRawFileName(),fMsqlDbHostStr(),
 fZebini(),fGcflag(),fRuninf(),flib_handle(),
 fJINFO(),fJFGES(),fJFDST(),fLQ(),fIQ(),fQ(),fdebug(),fRawFileExists()
{
  // method to be debugged and completed!
  CommInit();

  fdebug = debug;
  // Monte Carlo
  if(fdebug) cout << " --->debug TFndProcessRec::TFndProcessRec: "<< this->GetName() << " constructor called for MonteCarlo data" << endl;
  fTyp = new TString("MCAR");
  fNum = num;
  Char_t opt[10];
  sprintf(opt,"%d",num);
  SetFidaRC("NRUN",opt);
  if(fdebug) cout << " ---> ...opening library libFIDARC.so..." << endl; 
  flib_handle=dlopen("libFIDAMC.so",RTLD_LAZY);
  if(fdebug) cout << " --->done! Calling InitRun() " << endl; 
  
  InitRun(); 
  SetZebini();
  SetGcflag();
}

//_______________________________________________________________
TFndProcessRec::TFndProcessRec(const Char_t *typ, Int_t num, const Char_t *data_path, const Char_t *msql_db_host, Bool_t debug)
:TNamed("froot<->fidarc_interface","FINUDA froot<->rec_interface"),
 fProcDirBaseName(),fCurProcNum(),fTyp(),fNum(),fRawFileName(),fMsqlDbHostStr(),
 fZebini(),fGcflag(),fRuninf(),flib_handle(),
 fJINFO(),fJFGES(),fJFDST(),fLQ(),fIQ(),fQ(),fdebug(debug),fRawFileExists()
{ 

  CommInit();
  
  //--- Setting of $MSQL_DB_HOST
  if (msql_db_host && *msql_db_host == '$')
    msql_db_host = gSystem->Getenv(msql_db_host + 1);
  
  gSystem->Setenv("MSQL_DB_HOST",msql_db_host);
  if(fdebug) cout << "Starting \"FINUDA froot<->rec_interface\": msql host is: " << gSystem->Getenv("MSQL_DB_HOST") << endl;
  
  //--- Setting of $RDT
  if (data_path && *data_path == '$')
    data_path = gSystem->Getenv(data_path + 1);
  
  gSystem->Setenv("RDT",data_path);
  if(fdebug) cout << "Starting \"FINUDA froot<->rec_interface\": raw data path is: " << gSystem->Getenv("RDT") << endl;
  
  //--- Setting of temporary $INP (needed for running in multi-process mode)
  CreateProcInpFiles();
  
  // Real data
  //  fdebug = debug;
  fTyp = new TString(typ);
  fNum = num;
  fMsqlDbHostStr = new TString(msql_db_host);
  
  if(fdebug) cout << " --->debug TFndProcessRec::TFndProcessRec: "<< this->GetName() << " constructor called for Real data" << endl;
  CheckRawFile(); 
  
  TString type_pass;
  type_pass.Form("'%s'",fTyp->Data());
  
  if(fdebug) cout << " --->debug TFndProcessRec::TFndProcessRec: passing string " << type_pass << " to 'SetFidaRC' method" << endl;
  
  SetFidaRC("TRUN",(Char_t *)type_pass.Data());  
  
  Char_t opt[10];
  sprintf(opt,"%d",num);
  SetFidaRC("NRUN",opt);  
  
  // --- now chose the correct geometry according to the RUN-type/number
  //  ( really dangerous, but fortran needs it...the run-time from raw event should be used! )
  Bool_t IsNewDataTaking = kFALSE;
  Int_t RunNumThresholds[FROOT::RT_ONLM]; // last number of 2003 d-tak depending on the run type
  RunNumThresholds[FROOT::RT_FINU] = 2583;
  RunNumThresholds[FROOT::RT_COSM] = 4200; //check
  RunNumThresholds[FROOT::RT_PULS] = 1000;
  RunNumThresholds[FROOT::RT_PEDE] = 1000;
  RunNumThresholds[FROOT::RT_CALI] = 1000;
  RunNumThresholds[FROOT::RT_PLAT] = 1000;

  Int_t i_check;
  for(i_check=FROOT::RT_FINU; i_check<= FROOT::RT_ONLM; i_check++){
    //    Printf(" ---------------> \"%s\" <-> \"%s\"",fTyp->Data(),FROOT::RunType_Name(i_check).Data());
    if( !fTyp->CompareTo(FROOT::RunType_Name(i_check).Data()) ) break;
  }
  
  if(i_check!=FROOT::RT_ONLM){
    IsNewDataTaking = ( num <=  RunNumThresholds[i_check] ) ? kFALSE : kTRUE;
    Info("TFndProcessRec","Setting \"GEOC\" in \"fidarc.dat\" according to %d",IsNewDataTaking);
    (IsNewDataTaking) ? SetFidaRC("GEOC","1       1") : SetFidaRC("GEOC","0       1");
  }
  else{ // support for new data-taking only available in online mode
    Info("TFndProcessRec","run type is \"%s\": Setting \"GEOC\" in \"fidarc.dat\" to \"1\"");
    SetFidaRC("GEOC","1       1");
  }
  // ---
  SelectBhabhaInFidaRC(E_DOUBLE_HELIX); // by default use double-helix procedure for bhabha events

  // --- fidarc.dat card has been set

  flib_handle=dlopen("libFIDARC.so",RTLD_LAZY);

  if(fdebug) cout << "Initialyzing new run..." << endl;
  InitRun();
  if(fdebug) cout << "...done: now setting zebini..." << endl;
  SetZebini();
  if(fdebug) cout << "...done: now setting gcflag..." << endl;
  SetGcflag();
  if(fdebug) cout << "Constructor completed" << endl;
}

//_______________________________________________________________
TFndProcessRec::~TFndProcessRec()
{
  if(fdebug) cout << "TFndProcessRec debug: Destructor called." << endl; 
  
  // remove temporary (*.dat) directory
  //  gROOT->Info("~TFndProcessRec","DEBUG");
  CleanCurrentProcess();
  //  uglast_();
  //  gROOT->Info("~TFndProcessRec","DEBUG");
  delete fTyp;
  delete fMsqlDbHostStr;
  delete fRawFileName;
  //  gROOT->Info("~TFndProcessRec","DEBUG");
  fZebini = 0;
  fGcflag = 0;
  fRuninf = 0;
  flib_handle = 0;
  fIQ=0;
  fQ  =0;
  fLQ =0; 
  //  gROOT->Info("~TFndProcessRec","DEBUG");

  if(fdebug) cout << "TFndProcessRec debug: Destructor completed." << endl; 
}

//_______________________________________________________________
void TFndProcessRec::CommInit(){
  // Initialization needed by all ctors
  fProcDirBaseName = "fida_proc_";
}

//_______________________________________________________________
void TFndProcessRec::CreateProcInpFiles(){
  // make a directory for current process fida*.dat
  TString Ori_Inp = gSystem->Getenv("INP");
  // ---
  // Create a new subdirectory for current process setup files
  TString Proc_Inp = gSystem->Getenv("PWD");
  TString pth = Proc_Inp;

  Int_t i = 0;
  for(i=1; i <= K_NMAX_PROCESSES; i++){
    Proc_Inp.Form("%s/%s%d",pth.Data(),fProcDirBaseName.Data(),i);
    if(fdebug) Printf("looking for directory \"%s\"",Proc_Inp.Data());
    if( !gSystem->AccessPathName(Proc_Inp.Data()) ){
      Printf("Directory for process number \"%d\" already exists: trying next...",i);
      continue;
    }
    TString MkCmd = "mkdir ";
    MkCmd+=Proc_Inp;
    Printf("Sending command: \"%s\"",MkCmd.Data());
    gSystem->Exec(MkCmd.Data());
    break;
  }
  if(i>=K_NMAX_PROCESSES){
    Printf("Can not start process number \"%d\" (max # of processes reached)",i);
    FROOT::TerminateFroot(0,"TFndProcessRec -> Maximum number of processes reached");
  }
  fCurProcNum = i;
  
  // ---
  // copy every needed file in new directory
  TString copy_cmd;
  copy_cmd.Form("cp %s/*.dat %s/.",Ori_Inp.Data(),Proc_Inp.Data());
  Printf("Sending command: \"%s\"",copy_cmd.Data());
  gSystem->Exec(copy_cmd.Data());
  
  // --- 
  // change $INP for current session
  Printf("changing INP variable to directory \"%s\"",Proc_Inp.Data());
  gSystem->Setenv("INP",Proc_Inp.Data());
}

//_______________________________________________________________
void TFndProcessRec::CleanCurrentProcess(){
  // remove directory with current process fida*.dat

  TString Proc_Inp;
  TString Cdir = gSystem->Getenv("PWD");
  Proc_Inp.Form("%s/%s%d",Cdir.Data(),fProcDirBaseName.Data(),fCurProcNum);
  
  TString clean_cmd;
  clean_cmd.Form("rm -rf %s",Proc_Inp.Data());
  if(gSystem->AccessPathName(Proc_Inp.Data())){
    Warning("CleanCurrentProcess","Directory fpr current process does not exist (%s)",Proc_Inp.Data());
    return;
  }
  
  Printf( "Sending command: \"%s\"",clean_cmd.Data() );
  
  TString InpChk = gSystem->Getenv("INP");
  InpChk.Prepend("rm -rf ");
  Int_t CmpRes = InpChk.CompareTo(clean_cmd);
  if(fdebug) Printf("Comparing \"%s\" to \"%s\": result: %d",InpChk.Data(),clean_cmd.Data(),CmpRes);
  if(CmpRes != 0){
    Error("CleanCurrentProcess","$INP sems to be wrong!");
    FROOT::TerminateFroot();
  }
  gSystem->Exec(clean_cmd.Data());
}

//___________________________________
void TFndProcessRec::CheckRawFile()
{
  // while using TFndPreanMan this method should be
  // useless (raw-file check is performed previousely)

  if(fdebug) cout << "TFndProcessRec::CheckRawFile ---> fNum = " << fNum << endl;
  if(fNum==-1){// fNum==-1 means online (no raw file needed)
    fRawFileExists = kTRUE; // for any eventual check
    return;
  }

  fRawFileExists = kFALSE;
  // --- converting to 'run_num' to 'TString*'
  TString RunNumStr = "";
  RunNumStr.Form("%d",fNum);
  while(1){
    Int_t size = RunNumStr.Sizeof()-1;
    if(size>=5)break;
    RunNumStr.Prepend("0");
  }
  // --- filename and control strings
  delete fRawFileName;
  fRawFileName = new TString(fTyp->Data());
  *fRawFileName+=RunNumStr;
  if(fdebug) cout << endl << "Checking existence of run " << fRawFileName->Data();
  *fRawFileName+=".raw";
  // --- file must be in $RDT
  TString *PathName = new TString(gSystem->Getenv("RDT"));
  //cout << " in directory: " << PathName->Data() << endl;

  *PathName+="/";
  *PathName+=fRawFileName->Data();

  if(fdebug) cout << " --->debug TFndProcessRec::CheckRawFile(): file:" << PathName->Data() << endl; 
  // --- check raw file
  if(gSystem->AccessPathName(PathName->Data())==0){
   fRawFileExists = kTRUE;
   if(fdebug) cout << " ...raw file found"<< endl;
  }
  else{ // --- check zip-file
    *PathName+=".gz";
    if(gSystem->AccessPathName(PathName->Data())==0){
      fRawFileExists = kTRUE;
      if(fdebug) cout << " ...raw zip-file found"<< endl;
    }
  }
  // --- result
  if(!fRawFileExists) {
    cout << " *** file " << PathName->Data() << " not found!" << endl;
    cout << " Exiting program!" << endl << endl;
    gROOT->ProcessLine(".q");
  }
  delete PathName;
}

//_______________________________________________________________
void TFndProcessRec::InitRun() {
  // Method used to initialize fortran
  // see fidarc->main.f

  fMaxEv=100000000;

  Int_t a =  Ng;
  Int_t b = -Nh;  
  if(fdebug) cout << " --->debug TFndProcessRec::InitRun(): entering method." << endl; 

  gzcali_(&a);
  hlimit_(&b);
  uginit_();
    
  if(fdebug) cout << " --->debug TFndProcessRec::InitRun(): 'gzcali_', 'hlimit_' and 'uginit_' called succesfully." << endl; 
}  

//_________________________________________________________
void TFndProcessRec::SetZebini(){
  fZebini = (Zerbini_t *)dlsym(flib_handle,"zebini_");
  fLQ = &(fZebini->jinfo)-1;
  fIQ=fLQ+8;            // to get integer data
  fQ=(Float_t *)fIQ;    // to get integer data
  if(fdebug) cout << " ===>debug SetZebini(): fZebini address:" << fZebini << endl; 

}
//_________________________________________________________
void TFndProcessRec::SetGcflag(){
  fGcflag = (Gcflag_t *)dlsym(flib_handle,"gcflag_");
}


//_________________________________________________________
void TFndProcessRec::SetRuninf(){
  fRuninf = (Runinf_t *)dlsym(flib_handle,"runinf_");
}

//_______________________________________________________________
void TFndProcessRec::SelectBhabhaInFidaRC(E_prec_Bhabha_Mode mode){
  // uncomment or comment the BABA line according to the requested mode
  // mode: TFndProcessRec::E_prec_Bhabha_Mode { E_SINGLE_HELIX=0, E_DOUBLE_HELIX=1 };
  

  if(mode < E_SINGLE_HELIX || mode > E_DOUBLE_HELIX ){
    FROOT::TerminateFroot(0,"Mode not accepted in TFndProcessRec::SelectBhabhaInFidaRC");
  }
  fCurBhabhaRecMode = mode;

  if(fdebug) {
    cout << " --->debug TFndProcessRec::SelectBhabhaInFidaRC(): entering method" << endl;
    cout << "     mode  = " << mode << endl;
  }  
  TString fnam = FROOT::ExpandPathName("$INP"); 
  fnam+="/fidarc.dat";

  if(fdebug) cout << " --->debug TFndProcessRec::SetFidaRC(): reading file: " << fnam << endl;
  Char_t file[200][200];
  
  FILE *f = fopen(fnam.Data(),"r");
  int l=0;


  Int_t IsCardFormatOk = kFALSE;

  while (fgets(file[l],200,f)) {    Char_t *line=file[l++];
    line[strlen(line)-1]='\0';
    if(fdebug) cout << "=======> original line: \"" << line << "\"" << endl;
   
    if(strncmp("C BABA",line,6)==0){
      cout << "Current mode is \"double-helix\"" << endl;
      if(mode == E_SINGLE_HELIX){
        //      line[0]='\0';
        //      strcat(line,"BABA");
        strcpy(line,"BABA");
        cout << "fidarc card switched to single-helix mode" << endl;
      }
      IsCardFormatOk = kTRUE;
    }
    else if(strncmp("BABA",line,4)==0){
      cout << "Current mode is \"single-helix\"" << endl;
      if(mode == E_DOUBLE_HELIX){
        //      line[0]='\0';
        //      strcat(line,"C BABA");
        strcpy(line,"C BABA");
        cout << "fidarc card switched to double-helix mode" << endl;
      }
      IsCardFormatOk = kTRUE;
    }
    
    if(fdebug) cout << "=======> line after check: \"" << line << "\"" << endl;
  }
  fclose(f);
  
  if(IsCardFormatOk == kFALSE){
    Error("SelectBhabhaInFidaRC","(\"BABA\" or \"C BABA\" is expected somewhere) check spaces and endline!");
    FROOT::TerminateFroot(0,"Card format not recognized in TFndProcessRec::SelectBhabhaInFidaRC" );
  }
  // --- up to here file has not been modified yet  

  if(fdebug) cout << " --->debug TFndProcessRec::SetFidaRC(): modifying file: " << fnam << endl; 
  f = fopen(fnam.Data(),"w");
  for (int n=0;n<l;n++) fprintf(f,"%s\n",file[n]);
  fclose(f); 
  if(fdebug) cout << " --->debug TFndProcessRec::SetFidaRC(): done!" << endl; 
  
}

//________________________________________________________
void TFndProcessRec::SetFidaRC(Char_t *card, Char_t *option) {
  if(fdebug) {
    cout << " --->debug TFndProcessRec::SetFidaRC(): entering method" << endl;
    cout << "     card   = " << card << endl;
    cout << "     option = " << option << endl;
  }  
  Char_t fnam[200];
  strcpy(fnam,gSystem->Getenv("INP"));
  strcat(fnam,"/fidarc.dat");
  if(fdebug) cout << " --->debug TFndProcessRec::SetFidaRC(): reading file: " << fnam << endl;
  Char_t file[200][200];
  
  FILE *f = fopen(fnam,"r");
  int l=0;
  while (fgets(file[l],200,f)) {    Char_t *line=file[l++];
    line[strlen(line)-1]='\0';
    if (line[0]=='C') continue;
    
    int k=0;
    while (line[k]!=' ' && line[k]!='\0') k++;
    line[k]='\0';
    
    if (strcmp(card,line)==0) {
      strcat(line,"   ");
      strcat(line,option);
    }
    else line[k]=' ';
  }
  fclose(f);

  if(fdebug) cout << " --->debug TFndProcessRec::SetFidaRC(): modifying file: " << fnam << endl; 
  f = fopen(fnam,"w");
  for (int n=0;n<l;n++) fprintf(f,"%s\n",file[n]);
  fclose(f); 
  if(fdebug) cout << " --->debug TFndProcessRec::SetFidaRC(): done!" << endl; 

}

//_______________________________________________________________
void TFndProcessRec::InitNewRun(const Char_t *typ, Int_t num, const Char_t *data_path, const Char_t *msql_db_host, Bool_t debug){
  // start new run processing (recycling the used zebra structure)
  // (read_lib method called to open new raw file)

  //--- Setting of $MSQL_DB_HOST
  if(fdebug) cout << " --->debug TFndProcessRec::InitNewRun(): entering method." << endl; 

  if (msql_db_host && *msql_db_host == '$')
    msql_db_host = gSystem->Getenv(msql_db_host + 1);
  
  gSystem->Setenv("MSQL_DB_HOST",msql_db_host);
  if(fdebug) cout << "Starting \"FINUDA froot<->rec_interface\": msql host is: " << gSystem->Getenv("MSQL_DB_HOST") << endl;
  
  //--- Setting of $RDT
  if (data_path && *data_path == '$')
    data_path = gSystem->Getenv(data_path + 1);
  
  gSystem->Setenv("RDT",data_path);
  if(fdebug) cout << "Starting \"FINUDA froot<->rec_interface\": raw data path is: " << gSystem->Getenv("RDT") << endl;
  
  // Real data
  fdebug = debug;
  fTyp = new TString(typ);
  fNum = num;
  fMsqlDbHostStr = new TString(msql_db_host);
  
  if(fdebug) cout << " --->debug TFndProcessRec::TFndProcessRec: "<< this->GetName() << " constructor called for Real data" << endl;
  CheckRawFile(); 
  
  TString *type_pass = new TString(fTyp->Data());
  type_pass->Prepend("'");
  type_pass->Append("'");
  
  if(fdebug) cout << " --->debug TFndProcessRec::TFndProcessRec: passing string " << (Char_t *)type_pass->Data() << " to 'SetFidaRC' method" << endl;
  
  SetFidaRC("TRUN",(Char_t *)type_pass->Data());  
  delete type_pass;
  Char_t opt[10];
  sprintf(opt,"%d",num);
  SetFidaRC("NRUN",opt);  
  
  // ---
  fMaxEv=100000000;

  Int_t ierr = 0;
  Int_t IOUT[9];
  if(fNum!=-1){
    //cout << "===============> " << fRawFileName << endl;
    fRawFileName->Resize(9);
    //  TString finopen_str = TString(fRawFileName->Data());
    TString finopen_str = FROOT::ExpandPathName("$RDT");
    finopen_str+= "/";
    finopen_str+= fRawFileName->Data();
    //  cout << "DEBUG: Initializing new run with: passing string \"" << finopen_str.Data() << "\"" << endl;
    fin_open_((char *)finopen_str.Data(), IOUT, &ierr);
    //    fin_open03_((char *)finopen_str.Data(), IOUT, &ierr);
  }
  else fin_open_("online", IOUT, &ierr);
  
  if(!fRuninf) SetRuninf(); // check it out
  for(Int_t i=0;i<9;i++) fRuninf->iout[i] = IOUT[i];

}

//_______________________________________________________________
Int_t TFndProcessRec::UgEvent(){ 
  // Complete offline event processing (FORTRAN code called)
  // this method returns:
  //                     0 -> good event
  //                     1 -> run finished

  if(fdebug) cout << " --->debug TFndProcessRec::UgEvent(): entering method" << endl; 
  Int_t a=fZebini->ixdiv2;
  if(fdebug) cout << "----->UgEvent DEBUG1" << endl;  
  mzwipe_(&a);
  if(fdebug) cout << "----->UgEvent DEBUG2" << endl;  
  gtrigc_();
  if(fdebug) cout << "----->UgEvent DEBUG3" << endl;  
  gtrigi_();
  if(fdebug) cout << "----->UgEvent DEBUG4" << endl;  
  gutrev_();
  if(fdebug) cout << "----->UgEvent DEBUG5" << endl;
  
  if(fdebug) cout << " --->debug TFndProcessRec::UgEvent(): fortran called succesfully!" << endl;
  // zebra dst address assignment
  fJFDST=fZebini->jfdst;
  if(fdebug) cout << " ===>debug TFndProcessRec::fJFDST address: " << fJFDST << endl;
  
  
  if(fJFDST==0) {  // if 'fJFDST==0' run is finished or something is wrong
    if(fdebug)  cout << " ---> debug TFndProcessRec::UgEvent(): Run completed!" << endl << endl;
    return 1; // exit if run is finished
  }
  if(fdebug) cout << " --->debug TFndProcessRec::UgEvent(): event " <<  fIQ[fJFDST+2] << " procesed." << endl; 
  return 0;   
}

//_______________________________________________________________
UInt_t *TFndProcessRec::GetRawEvPointer(){
  return (UInt_t *) ( get_current_raw_event_() );
}

//_______________________________________________________________
/*
Int_t TFndProcessRec::UgEvent(){
  // this method returns:
  //                     0 -> good event
  //                     1 -> run completed

  if(FillJFGES()) return 1; // exit if run completed
  if(Reconstruct()) return 1;

  return 0;
}
*/

//_______________________________________________________________
Int_t TFndProcessRec::FillJFGES(){
  // Fill JFGES and evaluate detectors coordinates
  // without reconstructing the event
  // the event number (in fortran)
  // is automatically incremented by this method
  // this method returns:
  //                     0 -> good event
  //                     1 -> run completed
  //
  // see fidarc->gutrev

  if(fdebug) cout << " --->debug TFndProcessRec::FillJFGES(): entering method" << endl;
  //  cout << "DEBUG_FILLGES 1: Runinf address: " << fRuninf << endl;

  Int_t a=fZebini->ixdiv2;
  //  cout << "DEBUG_FILLGES 2: a = " << a << " --- fZebini = " << fZebini << endl;
  mzwipe_(&a);
  gtrigc_();
  gtrigi_();

  // --- GUTREV first part
  Int_t ierr = 0;
  rdtupk_(&ierr); // fill zebra GES
  if(ierr < 0) return 1; // run completed

  checktrig_(); // added into this code from fidarc-version 601 on

  //  consinit_(); // (from v521 on) initialize constant in fidarc
  //    consinit has been included into uginit from fidarc-version 601 on

  vrecon_();  // reconstruct vertex
  films_();   // needed by microstrips
  patstrw_(); // evaluate straw-tube triplets
  //  provatof_(); // (used by TOF group for dedicated analysis) 
  // ---

  fJINFO=fZebini->jinfo; // zebra info addresses assignment (event-zero information)
  fJFGES=fZebini->jfges; // zebra ges addresses assignment
  
  return 0;
  //return fGcflag->ieorun;
} 

//_______________________________________________________________
Int_t TFndProcessRec::Reconstruct(){
  // reconstruct current event (FORTRAN code)
  // This method corresponds to the second part of GUTREV (fidarc fortran code)
  // NOTE: This method must be called only if FillJFGES()
  // was called previousely
  // this method returns:
  //                     0 -> good event
  //                     1 -> run completed
  //                    -1 -> problems
  //  
  //  COSM NOT IMPLEMENTED
  //  see fidarc->gutrev_() from "IF (HYPE) THEN" on [without rdtupk]

  fJFDST = 0;
  if(!fJFGES){
    Warning("Reconstruct","fJFGES not assigned.");
    return -1;
  }
  
  Int_t trig_type = fIQ[fJFGES+3];
  if(fdebug) cout << "trigger type: " << trig_type << endl;

  // the following code is the "translation" of the GUTREV fortran routine
  switch(trig_type){
  case 1: // HYPE
    genrec_();
    trkfit_();
#ifdef FIDA_USE_SHORT
    ktrkvert_(); // from fida-version 603 on
    vertfit_();
#endif
#ifdef FIDA_USE_KALMAN
    kalmanfit_();
#endif
    neutfind_();
    break;
  case 2: // BHABHA
    if(fCurBhabhaRecMode == E_DOUBLE_HELIX){
      bhabharc_();
      bhatrkfit_();
    }
    else if(fCurBhabhaRecMode == E_SINGLE_HELIX){
      babarc_();
      babafit_();
    }
    break;
  case 3: // COSM
    Warning("Reconstruct()","COSM reconstruction not implemented");
    return 0;
  case 11: // HYPE && BHABHA
    genrec_();
    trkfit_();
#ifdef FIDA_USE_SHORT
    ktrkvert_(); // from fida-version 603 on
    vertfit_();
#endif
#ifdef FIDA_USE_KALMAN
    kalmanfit_();
#endif
    neutfind_();
    break;
  case 21:
    Warning("Reconstruct()","MONTECARLO reconstruction not implemented");
    return 0;
  default: return 0;

  }

  fJFDST=fZebini->jfdst;  // zebra dst addresses assignment

  if(fdebug) cout << " ===>debug TFndProcessRec::fJFDST address: " << fJFDST << endl;
  if(fJFDST==0) {  // if 'fJFDST==0' run is finished or something is wrong
    if(fdebug)  cout << " ---> debug TFndProcessRec::Reconstruct(): Run completed!" << endl << endl;
    return 1; // exit if run is finished
  }
  if(fdebug) cout << " --->debug TFndProcessRec::Reconstruct(): event " <<  fIQ[fJFDST+2] << " procesed." << endl; 
  return 0;   

} 

//_______________________________________________________________
void TFndProcessRec::CloseRun(Bool_t last) {
  // Close raw file (read_lib method called)
  // Perform FORTRAN close operation (uglast/zebra reset) if last is true
  // Should be called with "last=true" before calling CloseInterface()

  if(fdebug) cout << "TFndProcessRec::CloseRun(); Closing raw file..." << endl;
  if(!last){
    if(fdebug) cout << "--> TFndProcessRec::CloseRun() DEBUG 1" << endl;
    Int_t r_tim = fIQ[fJFGES+4];
    if(!r_tim) r_tim = fIQ[fJFDST+4];
    if(fdebug) cout << "Run-time is: " << r_tim <<  endl;
    if(!r_tim) return;
    if(r_tim <= FROOT::DTak_2006_SftwStartTime) fin_close03_(); 
    else fin_close_();

    if(fdebug) cout << "TFndProcessRec::CloseRun(); raw file closed: Runinf address: " << fRuninf << endl;
  }
  else{
    if(fdebug) cout << "--> TFndProcessRec::CloseRun() DEBUG 2" << endl;
    uglast_();
    if(fdebug) cout << "--> TFndProcessRec::CloseRun() DEBUG 3" << endl;
    Int_t a;
    a = fZebini->ixdiv1;   mzwipe_(&a);
    if(fdebug) cout << "--> TFndProcessRec::CloseRun() DEBUG 4" << endl;
    a = fZebini->ixdiv2;   mzwipe_(&a);
    if(fdebug) cout << "--> TFndProcessRec::CloseRun() DEBUG 5" << endl;
    a = fZebini->ixdiv3;   mzwipe_(&a);
    if(fdebug) cout << "--> TFndProcessRec::CloseRun() DEBUG 6" << endl;
  }
  CleanRecOutputs(kTRUE,kTRUE,kTRUE);
}

//___________________________________
Int_t TFndProcessRec::CloseInterface(){
  // close fida*.so
  // Should be called before calling dtor
  // return value: exit code of the "dlclose" function


  Int_t dlclose_err = dlclose(flib_handle);
  if(fdebug) cout << "dlclose output: " << dlclose_err << endl;
  return dlclose_err;
}

//_________________________________________________________
Bool_t TFndProcessRec::GetStatusBit(Int_t l, Int_t k) {
  // this function returns the value of the bit 'k' of the bank 'l'
  // (grab bit number "k" from the status word of bank "l")
  // range for "k": 1-2-3-4-5-...

  if(k<=0){
    Error("GetStatusBit","bit number must be >0 (%d received)",k);
    FROOT::TerminateFroot();
    return 0;
  }

  if (fJFDST==0) return 0;
  else return (Bool_t)(fIQ[l]&((Int_t)TMath::Power(2,k-1)));
}

//___________________________________________________________
void TFndProcessRec::PrintZebiniContent(){
  if(fZebini){
    cout << "fZebini->ixdst: " << fZebini->ixdst << "   at address: " << &(fZebini->ixdst) << endl;
    cout << "fZebini->ixdiv1: " << fZebini->ixdiv1 << "   at address: " << &(fZebini->ixdiv1) << endl;
    cout << "fZebini->ixdiv2: " << fZebini->ixdiv2 << "   at address: " << &(fZebini->ixdiv2) << endl;
    cout << "fZebini->ixdiv3: " << fZebini->ixdiv3 << "   at address: " << &(fZebini->ixdiv3) << endl;
    for(Int_t i=0;i<16;i++)
      cout << "fZebini->fence[" <<i << "]: " << fZebini->fence[i] << "   at address: " << &(fZebini->fence[i]) << endl;
    cout << "fZebini->jinfo: " << fZebini->jinfo << "   at address: " << &(fZebini->jinfo) << endl;
    cout << "fZebini->jfges: " << fZebini->jfges << "   at address: " << &(fZebini->jfges) << endl;
    cout << "fZebini->jmstr: " << fZebini->jmstr << "   at address: " << &(fZebini->jmstr) << endl;
    cout << "fZebini->jstrw: " << fZebini->jstrw << "   at address: " << &(fZebini->jstrw) << endl;
    cout << "fZebini->jfdst: " << fZebini->jfdst << "   at address: " << &(fZebini->jfdst) << endl;
  }
  else cout << "No zebini defined!!!" << endl;
}

//___________________________________________________________
void TFndProcessRec::PrintGcflagContent(){
  if(fGcflag){
    cout << "fGcflag->idebug ===> " << fGcflag->idebug << "   at address: " << &(fGcflag->idebug) << endl;
    cout << "fGcflag->idemin ===> " << fGcflag->idemin << "   at address: " << &(fGcflag->idemin) << endl;
    cout << "fGcflag->idemax ===> " << fGcflag->idemax << "   at address: " << &(fGcflag->idemax) << endl;
    cout << "fGcflag->itest ====> " << fGcflag->itest << "   at address: " << &(fGcflag->itest) << endl;
    cout << "fGcflag->idrun ====> " << fGcflag->idrun << "   at address: " << &(fGcflag->idrun) << endl;
    cout << "fGcflag->idevt ====> " << fGcflag->idevt << "   at address: " << &(fGcflag->idevt) << endl;
    cout << "fGcflag->ieorun ===> " << fGcflag->ieorun << "   at address: " << &(fGcflag->ieorun) << endl;
    cout << "fGcflag->ieotri ===> " << fGcflag->ieotri << "   at address: " << &(fGcflag->ieotri) << endl;
    cout << "fGcflag->ievent ===> " << fGcflag->ievent << "   at address: " << &(fGcflag->ievent) << endl;
    for(Int_t i=0;i<10;i++)
      cout << "fGcflag->iswit " << i << " ==> " << fGcflag->iswit[i] << "   at address: " << &(fGcflag->iswit[i]) << endl;
    for(Int_t i=0;i<20;i++)
      cout << "fGcflag->ifinit " << i << " => " << fGcflag->ifinit[i] << "   at address: " << &(fGcflag->ifinit[i]) << endl;
    cout << "fGcflag->nevent ===> " << fGcflag->nevent << "   at address: " << &(fGcflag->nevent) << endl;
    for(Int_t i=0;i<2;i++)
      cout << "fGcflag->nrndm " << i << " ==> " << fGcflag->nrndm[i] << "   at address: " << &(fGcflag->nrndm[i]) << endl;
  }
  else cout << "No gcflag defined!!!" << endl;
}

//___________________________________________________________
void TFndProcessRec::CleanRecOutputs(Bool_t out,Bool_t his,Bool_t dst){

  Printf("TFndProcessRec::CleanRecOutputs ---> Cleaning reconstruction outputs.");
  TString dest_cmd_out;
  TString dest_cmd_his;
  TString dest_cmd_dst;
  
  dest_cmd_out.Form("rm -f %s/*.lis", FROOT::ExpandPathName("$REC_OUTPUTS/out").Data() );
  dest_cmd_his.Form("rm -f %s/*.rzf", FROOT::ExpandPathName("$REC_OUTPUTS/histos").Data() );
  dest_cmd_dst.Form("rm -f %s/*.dat", FROOT::ExpandPathName("$REC_OUTPUTS/dstapes").Data() );

  if(out){
    Printf(" ---> Sending command \"%s\"",dest_cmd_out.Data());
    gSystem->Exec(dest_cmd_out.Data());
  }
  if(his){
    Printf(" ---> Sending command \"%s\"",dest_cmd_his.Data());
    gSystem->Exec(dest_cmd_his.Data());
  }
  if(dst){
    Printf(" ---> Sending command \"%s\"",dest_cmd_dst.Data());
    gSystem->Exec(dest_cmd_dst.Data());
  }


  Printf("TFndProcessRec::CleanRecOutputs ---> Reconstruction outputs cleaned.");
}

// ---------- EXAMPLE methods -----------
//___________________________________________________________
void TFndProcessRec::AutoProcessRun() {
  // This method must not be used in froot: is just an example
  fBhaFlag  = kFALSE;
  fHypeFlag = kFALSE;
  fCosmFlag = kFALSE;

  Int_t nev=0;
  if(fdebug) cout << " --->debug ProcessRun(): beginning events loop." << endl; 
  
  Int_t end_of_run=0;
  while(end_of_run==0 && nev<fMaxEv){
    if(nev==fMaxEv) cout << "Max number of events reached:" << endl;
    nev++;
    cout << "Start processing event " << nev << endl;
    Int_t a=fZebini->ixdiv2;
    mzwipe_(&a);
    gtrigc_();
    gtrigi_();
    gutrev_();
    //---------------------------
    fJFDST=fZebini->jfdst;
 
    if(fdebug) cout << " ===>debug ProcessRun():fJFDST = " << fJFDST << endl;
    if (fJFDST==0) {
      cout << "fJFDST not assigned to anything!!! exiting from loop" << endl;
      end_of_run = 1;
    }
    // here I begin to read zebra dst-bank:
    // reading run number and event number (using fJFDST)
    cout << "run number: " << fIQ[fJFDST+1]<< endl;
    cout << "ev. number: " << fIQ[fJFDST+2]<< endl;
    Int_t trig_type = fIQ[fJFDST+3];
    cout << "trigger type: " << trig_type << endl;
    /*
      Values of trigger bit:
       1 ---> Hypernuclear trigger
       2 ---> Bhabha Trigger
       3 ---> Cosmic ray trigger
      11 ---> Both HYP and BHA bits on (priority to HYP trigger)
      21 ---> Single track (only MC)
    */
    switch (trig_type) {
    case 1:
      cout << " -> trigger: HYP" << endl;
      end_of_run = GetHypeDst();
      break;
    case 2:
      cout << " -> trigger: BHA" << endl;
      end_of_run = GetBabaDst();
      //end_of_run = GetBhabhaDst();
      break;
    case 3:
      cout << " -> trigger: COSM" << endl;
      end_of_run = GetCosmDst();
      break;
    case 11:
      cout << " -> double trigger: HYP & BHA (reading HYPE bank)" << endl;
      end_of_run = GetHypeDst();
      break;
    case 21:
      cout << " -> Single track (Montecarlo only)" << endl;
      end_of_run = 1;
      break;
    default: 
      cout << "unknown trigger bit value!!!" << endl;
    }

    if(end_of_run) cout << " This event doesn't exist: run is finished." << endl;
  }
  if(fdebug) cout << " --->debug ProcessRun(): events loop completed" << endl; 
  cout << "Processing completed." << endl; 

  // Drawing histos...
  TFile *OutRootFile = new TFile("Autoprocess_histos.root","RECREATE");
  //--
  if(fHypeFlag){
    TCanvas *hype_canv = new TCanvas("hype_canv","hype canvas",20,20,800,800);
    hype_canv->Divide(2,2);
    hype_canv->cd(1);
    fHyp_K_min_xy_vert_histo->Draw();
    fHyp_K_min_xy_vert_histo->Write();
    hype_canv->cd(2);
    fHyp_K_min_z_vert_histo->Draw();
    fHyp_K_min_z_vert_histo->Write();
    hype_canv->cd(3);
    fHyp_K_min_mom_neg_trk->Draw();
    fHyp_K_min_mom_neg_trk->Write();
    hype_canv->cd(4);
    fHyp_K_min_mom_pos_trk->Draw();
    fHyp_K_min_mom_pos_trk->Write();

    cout << fcount_k_min << " Hypernuclear events ; " << fcount_k_min_inside << " stopped in target (" << 100 * ((Float_t)fcount_k_min_inside / (Float_t)fcount_k_min) << "%)" << endl;
  }
  //--
  if(fBhaFlag){
    TCanvas *bhabha_canv = new TCanvas("bhabha_canv","bhabha canvas",40,40,1000,500);
    bhabha_canv->Divide(2,1);
    bhabha_canv->cd(1);
    fBABA_mom_patt_histo->Draw();
    bhabha_canv->cd(2);
    fBABA_mom_fit_histo->Draw();
    cout << fcount_BABA << " Bhabha events ; " << fcount_BABA_OK << " fitted events (" << 100 * ((Float_t)fcount_BABA_OK / (Float_t)fcount_BABA) << "%)" << endl;
  }
  //--
  if(fCosmFlag){
    TCanvas *cosm_canv = new TCanvas("cosm_canv","cosm canvas",60,60,900,300);
    cosm_canv->Divide(3,1);
    cosm_canv->cd(1);
    fCOSM_mom_patt_histo->Draw();
    cosm_canv->cd(2);
    fCOSM_mom_fit_histo->Draw();
    cosm_canv->cd(3);
    fCosm_SignOfCharge_histo->Draw();
   
    cout << fcount_COSM << " Cosmic events ; " << fcount_COSM_OK << " fitted events (" << 100 * ((Float_t)fcount_COSM_OK / (Float_t)fcount_COSM) << "%)" << endl;
  }
  //--
  OutRootFile->Close();
}

//_____________________________
Int_t TFndProcessRec::GetHypeDst() {
  // This method must not be used in froot: is just an example
  // ---> k- bank access
  if(fJFDST==0) { cout << "no fJFDST!!!" << endl; return 1; }
  
  if(fHypeFlag==kFALSE){ // Declare histos only the first time
    fHyp_K_min_xy_vert_histo = new TH2F("K_min_xy_vert_histo","K- vertex: x-y coordinates",200,-10,10,200,-10,10);
    fHyp_K_min_z_vert_histo = new TH1F("K_min_z_vert_histo","K- vertex: z coordinate",200,-10,10);
    fHyp_K_min_mom_neg_trk = new TH1F("fHyp_K_min_mom_neg_trk","K- negative trk: momentum at vertex",1400,0,0.7);
    fHyp_K_min_mom_pos_trk = new TH1F("fHyp_K_min_mom_pos_trk","K- positive trk: momentum at vertex",1000,0,2);
    fcount_k_min = 0;
    fcount_k_min_inside = 0;
    fHypeFlag=kTRUE;
  }
  fcount_k_min++;
  Int_t LKMIN = fLQ[fJFDST-1];
  Int_t Interact_inside = GetStatusBit(LKMIN,1); // 1==>Interaction inside the target
  if (LKMIN)  {         
    //cout << "LKMIN offset= " << LKMIN << endl;
    cout << "vertex ======> " << fQ[LKMIN+5] << " - " << fQ[LKMIN+6] << " - "<< fQ[LKMIN+7] << endl;
    fHyp_K_min_xy_vert_histo->Fill(fQ[LKMIN+5],fQ[LKMIN+6]);  
    fHyp_K_min_z_vert_histo->Fill(fQ[LKMIN+7]); 
    // K- negative track bank access:
    Int_t LKMIN_NEGT = fLQ[LKMIN-3];
    if(LKMIN_NEGT){
      Int_t fit_err_code = fIQ[LKMIN_NEGT+7];
      if(fit_err_code==0 && Interact_inside==1){
        fHyp_K_min_mom_neg_trk->Fill(fQ[LKMIN_NEGT+8]);
        fcount_k_min_inside++;
      }
    }
    // K- positive track bank access:
    Int_t LKMIN_POST = fLQ[LKMIN-4];
    if(LKMIN_POST){
      Int_t fit_err_code = fIQ[LKMIN_POST+7];
      if(fit_err_code==0 && Interact_inside==1)
        fHyp_K_min_mom_pos_trk->Fill(fQ[LKMIN_POST+8]);  
    }
  }
  cout << " - - - - - - - " << endl << endl;
  return 0; // return '0' if event is exists
}

//_____________________________
Int_t TFndProcessRec::GetBabaDst() {
  // This method must not be used in froot: is just an example
  // baba bank access
  // Shows reconstructed momentum for Bhabha events (simplified routine)
  
  if(fJFDST==0) { cout << "no fJFDST!!!" << endl; return 1; }
  
  if(fBhaFlag==kFALSE){ // Declare histos only the first time
    fBABA_mom_patt_histo = new TH1F("fBABA_mom_patt_histo","BABA bank: momentum from pattern recognition",1000,0,1);
    fBABA_mom_fit_histo = new TH1F("fBABA_mom_fit_histo","BABA bank: momentum from fit",1000,0,1);
    fcount_BABA=0;
    fcount_BABA_OK=0;
    fBhaFlag=kTRUE;
  }
  fcount_BABA++;
  Int_t Fbaba = fLQ[fJFDST-1];
  if(Fbaba){
    Int_t FitOk = GetStatusBit(Fbaba,1);
    if(FitOk){
      fBABA_mom_patt_histo->Fill(fQ[Fbaba+7]);
      fBABA_mom_fit_histo->Fill(fQ[Fbaba+10]);
      fcount_BABA_OK++;
    }    
  }   
  
  cout << " - - - - - - - " << endl << endl;
  return 0; // return '0' if event is exists
}

//_____________________________
Int_t TFndProcessRec::GetBhabhaDst() {
  // This method must not be used in froot: is just an example
  // bhabha bank access

  if(fJFDST==0) { cout << "no fJFDST!!!" << endl; return 1; }


  cout << " - - - - - - - " << endl << endl;
  return 0; // return '0' if event is exists
}

//_____________________________
Int_t TFndProcessRec::GetCosmDst() { // B>0
  // This method must not be used in froot: is just an example
  if(fJFDST==0)  { cout << "no fJFDST!!!" << endl; return 1; }
  
  if(fCosmFlag==kFALSE){ // Declare histos only the first time
    fCosm_SignOfCharge_histo = new TH1F("fCosm_SignOfCharge_histo","Cosmic DST: Sign of charge Histogram",3,-1,1);
    fCOSM_mom_patt_histo = new TH1F("fCOSM_mom_patt_histo","COSM bank: momentum from pattern recognition",100,0,6);
    fCOSM_mom_fit_histo = new TH1F("fCOSM_mom_fit_histo","COSM bank: momentum from fit",100,0,6);
    fcount_COSM=0;
    fcount_COSM_OK=0;
    fCosmFlag=kTRUE;
  }
  fcount_COSM++;
  Int_t Fcosm = fLQ[fJFDST-1];
  if(Fcosm){
    Int_t Fit_OK = GetStatusBit(Fcosm,1);
    if(Fit_OK==0){
      fCosm_SignOfCharge_histo->Fill(fIQ[Fcosm+1]);
      fCOSM_mom_patt_histo->Fill(fQ[Fcosm+10]);
      fCOSM_mom_fit_histo->Fill(fQ[Fcosm+11]);
      fcount_COSM_OK++;      
    }
  }
  cout << " - - - - - - - " << endl << endl;
  return 0; // return '0' if event is exists
}

---