CalorimeterLevel2/src/CaloCore.cpp

//
//   Given a calibration and a data file this program create an ntuple with LEVEL2 calorimeter variables - Emiliano Mocchiutti
//
//   CaloCore.cxx      version 3.05 (2006-05-30)
//
//   The only input needed is the path to the directory created by YODA for the data file you want to analyze.
//
//   Changelog:
//
//   3.08 (2006-11-13): Added high energy nuclei capability and "process all events" capability.
//
//   3.04 - 3.05 (2006-05-30): Qlast and nlast are now calculated using 4 (not 8) strips aournd the shower axis. Small bug fixed.
//
//   3.03 - 3.04 (2006-05-23): Forgot to put impx and impy in the PAMELA reference system, fixed.
//
//   3.02 - 3.03 (2006-05-18): updated to be called in DarthVader. Output dimension are now in cm and in the PAMELA reference system.
//
//   3.01 - 3.02 (2006-04-21): when copying entries get size of caclone and not of ca... (shouldn't matter). Fixed increasing file dimension bug when reprocessing.
//                             Added variable planemax[2], plane of maximum energy release (x and y) in final output. Use ItoRunInfo instead of RunInfo.
//
//   3.00 - 3.01 (2006-04-14): fixed small bug in tagging the track used to determine track-related variables, put in caloprocessing the opening of parameters files, fixed
//                             small bug in fortran routines
//
//   2.01 - 3.00 (2006-04-14): almost everything has changed. Now it can process one, all or some runs, introduced the final CaloLevel2 class+methods and the 
//                             working class "CaloProcessing", linked to the preliminary tracker flight software v0r00, reads YODA unique output files,
//                             reduced the number of installed libraries, F77 programs splitted depending on the function contained, introduced the readout and
//                             processing of self-trigger events, if the tracker provides more than one track all calorimeter track-related variables are saved
//                             as many times as the number of tracks (via the TClonesArray object in the level2 rootple) and many other small changes.
//
//   2.00 - 2.01 (2006-01-26): bug: wrong calculation of baselines in some cases, fixed.
//
//   1.00 - 2.00 (2006-01-11): use TSQL ROOT classes instead of directly calling MySQL.
//
//   0.00 - 1.00 (2005-09-14): seems working.
//
//   0.00        (2005-09-09): clone of CaloLEVEL2.c .
//
// C/C++ headers
//
#include <fstream>
#include <string.h>
//
// ROOT headers
//
#include <TTree.h>
#include <TClassEdit.h>
#include <TObject.h>
#include <TList.h>
#include <TArrayI.h>
#include <TSystem.h>
#include <TSystemDirectory.h>
#include <TString.h>
#include <TFile.h>
#include <TClass.h>
#include <TCanvas.h>
#include <TH1.h>
#include <TH1F.h>
#include <TH2D.h>
#include <TLatex.h>
#include <TPad.h>
#include <TSQLServer.h>
#include <TSQLRow.h>
#include <TSQLResult.h>
#include <TClonesArray.h>
#include <TStreamerInfo.h>
//
// YODA headers
//
#include <PamelaRun.h>
#include <physics/trigger/TriggerEvent.h>
//
// This program headers
//
#include <RunInfo.h>
#include <CaloCore.h>
#include <CaloLevel2.h>
#include <CaloProcessing.h>
#include <CaloVerl2.h>
//
// Tracker classes headers and definitions
//
#include <TrkLevel2.h>
//
using namespace std;
//
// CORE ROUTINE
//
int CaloCore(UInt_t run, TFile *file, TSQLServer *dbc, Int_t calargc, char *calargv[]){
  //
  // Set these to true to have a very verbose output.
  //
  Bool_t verbose = false;
  Bool_t debug = false;
  //
  Bool_t trackanyway = true;
  //
  Float_t rigdefault = 50.;
  //
  Bool_t hZn = true;
  //
  Bool_t withtrk = true;
  //
  Bool_t st = true;
  //
  // Output directory is the working directoy.
  //
  const char* outdir = gSystem->DirName(gSystem->DirName(file->GetPath()));
  //
  Int_t ri = 0;
  TString processFolder = "calorimeterFolder";
  if ( calargc > 0 ){
    ri = 0;
    while ( ri < calargc ){
      if ( !strcmp(calargv[ri],"-processFolder") ) {
        if ( calargc < ri+1 ){
          throw -3;
        };
        processFolder = (TString)calargv[ri+1]; 
        ri++;
      };
      if ( !strcmp(calargv[ri],"-v") ||  !strcmp(calargv[ri],"--verbose") ) {
        verbose = true;
      };
      if ( !strcmp(calargv[ri],"-g") ||  !strcmp(calargv[ri],"--debug") ) {
        debug = true;
      };
      if ( !strcmp(calargv[ri],"--alltracks")  ) {
        trackanyway = true;
      };
      if ( !strcmp(calargv[ri],"--no-tracker")  ) {
        withtrk = false;
      };
      if ( !strcmp(calargv[ri],"--with-tracker")  ) {
        withtrk = true;
      };
      if ( !strcmp(calargv[ri],"--defrig")  ) {
        if ( calargc < ri+1 ){
          throw -3;
        };
        rigdefault = atof(calargv[ri+1]);       
        ri++;
      };
      if ( !strcmp(calargv[ri],"--no-alltracks")  ) {
        trackanyway = false;
      };
      if ( !strcmp(calargv[ri],"--highZnuclei")  ) {
        hZn = true;
      };
      if ( !strcmp(calargv[ri],"--no-highZnuclei")  ) {
        hZn = false;
      };
      if ( !strcmp(calargv[ri],"--selftrigger")  ) {
        st = true;
      };
      if ( !strcmp(calargv[ri],"--no-selftrigger")  ) {
        st = false;
      };
      if ( !strcmp(calargv[ri],"--help")  ) {
        printf("\n\n CALORIMETER HELP CALLED\n\n");
        printf(" CaloCore options: \n");
        printf("  -v | --verbose       be verbose\n");
        printf("  -g | --debug         be really verbose\n");
        printf("  --defrig rig         rig is the default rigidity in GV to be used to\n");
        printf("                       obtain calorimeter variables in the routines\n");
        printf("                       \"alltracks\" and \"higZnuclei\" [default = 50]\n");
        printf("  --alltracks          fill the track related variables even in the case\n");
        printf("                       of no tracks from tracker and no selftrigger event\n");
        printf("                       when we have a calorimeter fit for both views\n");
        printf("  --no-alltracks        fill the track related variables only in the case\n");
        printf("                       of a good track from tracker or selftrigger [default]\n");
        printf("  --highZnuclei        call the routine to analyze high Z nuclei\n");
        printf("                       selftrigger events [default]\n");
        printf("  --no-highZnuclei     do not call the routine to analyze high Z nuclei\n");
        printf("                       selftrigger events\n");
        printf("  --no-tracker         do not use tracker level2\n"); 
        printf("  --with-tracker       use tracker level2 [default]\n"); 
        printf("  --selftrigger        process selftrigger events [default]\n");
        printf("  --no-selftrigger     skip selftrigger events\n");
        throw -114;
      };
      ri++;
    };
  };
  //
  if ( verbose ) printf("\n");
  if ( hZn && verbose ) printf(" Calling high energy nuclei subroutine \n");
  if ( trackanyway && verbose ) printf(" Filling track related variables for all the possible tracks \n");
  if ( verbose && ( hZn || trackanyway ) ) printf(" Default assumed rigidity %f \n",rigdefault);
  if ( verbose && st ) printf(" Calling selftrigger subroutine \n");
  if ( verbose && withtrk ) printf(" Using tracker level2 data \n");
  //
  // Working filename
  //
  TString outputfile;
  stringstream name;
  name.str("");
  name << outdir << "/";
  //
  // Variables.
  //
  //  TFile *file = 0;
  TTree *tracker = 0;
  TTree *calo = 0;
  TTree *caloclone = 0;
  Bool_t reproc = false;
  Bool_t reprocall = false;
  UInt_t nevents = 0;
  UInt_t nobefrun = 0;
  UInt_t noaftrun = 0;
  UInt_t numbofrun = 0;
  UInt_t totnorun = 0;
  //
  Int_t code = 0;
  Int_t sgnl;
  //
  // calorimeter level2 classes
  //
  CaloLevel2 *ca = new CaloLevel2();
  CaloLevel2 *caclone = new CaloLevel2();
  //
  TrkLevel2 *trk = new TrkLevel2();
  Int_t nevtrkl2 = 0;
  //
  UInt_t procev = 0;
  //
  // define variables where to store the absolute run header and run trailer times (unsigned long long integers, when set to a number use to store the correct number).
  //
  UInt_t runheadtime = 0;
  UInt_t runtrailtime = 0;
  UInt_t evfrom = 0;
  UInt_t evto = 0;
  UInt_t totfileentries = 0;
  UInt_t idRun = 0;
  Int_t id_reg_run=-1;
  stringstream ftmpname;
  TString fname;
  //
  // define variables for opening and reading level0 file
  //
  TFile *l0File = 0;
  TTree *l0tr = 0;
  TBranch *l0head = 0;
  TBranch *l0calo = 0;
  TBranch *l0trig = 0;
  pamela::EventHeader *eh = 0;
  pamela::PscuHeader *ph = 0;
  pamela::trigger::TriggerEvent *trig = 0;
  //
  // Define some basic variables
  // 
  CaloProcessing *event = new CaloProcessing(); // NOTICE: very important to call here the constructor!
  stringstream file2;
  stringstream file3;
  stringstream qy;
  //  Bool_t imtrack = false;
  Bool_t filled = false;
  //
  UInt_t caloevents = 0;
  stringstream calfile;
  stringstream aligfile;
  //
  Int_t i = -1;    
  Int_t itr = -1;    
  Int_t badevent = 0;
  Int_t totevent = 0;
  //
  UInt_t atime = 0;
  //
  Int_t S3 = 0;
  Int_t S2 = 0;
  Int_t S12 = 0;
  Int_t S11 = 0;
  UInt_t re = 0;
  UInt_t jumped = 0;
  //
  TString caloversion;
  ItoRunInfo *runinfo = 0;
  TArrayI *runlist = 0;
  //
  Float_t tmptrigty = -1.;      
  Int_t ntrkentry = 0;
  GL_PARAM *q4 = new GL_PARAM();
  UInt_t tttrkpar1 = 0;
  Bool_t trkpar1 = true;
  GL_ROOT *glroot = new GL_ROOT();
  GL_TIMESYNC *dbtime = 0;
  //
  //
  //
  TFile *tempfile = 0;
  TTree *tempcalo = 0;
  stringstream tempname;
  stringstream calofolder;
  tempname.str("");
  tempname << outdir;
  tempname << "/" << processFolder.Data();
  calofolder.str("");
  calofolder << tempname.str().c_str();
  gSystem->MakeDirectory(calofolder.str().c_str());
  tempname << "/calotree_run";
  tempname << run << ".root";  
  //
  // As a first thing we must check what we have to do: if run = -1 we must process all events in the file has been passed
  // if run != -1 we must process only that run but first we have to check if the branch calorimeter already exist in the file
  // if it exists we are reprocessing data and we must delete that entries, if not we must create it.
  // 
  if ( run == 0 ) reproc = true;
  //
  //
  //
  if ( !file->IsOpen() ){
    if ( verbose ) printf(" CALORIMETER - ERROR: cannot open file for writing\n");
    throw -101;    
  };
  //
  if ( withtrk ){
    //
    // Does it contain the Tracker tree?
    //
    tracker = (TTree*)file->Get("Tracker");
    if ( !tracker ) {
      if ( verbose ) printf(" CALORIMETER - ERROR: no tracker tree\n");
      code = -102;
      goto closeandexit;
    };
    tracker->SetBranchAddress("TrkLevel2",&trk);  
    nevtrkl2 = tracker->GetEntries();
  };
  //
  // Call runinfo
  //  
  sgnl = 0;
  runinfo = new ItoRunInfo(file);
  //
  // update versioning information and retrieve informations about the run to be processed
  //
  caloversion = CaloInfo(false);
  sgnl = runinfo->Update(run,"CALO",caloversion);
  //
  if ( sgnl ){
    if ( verbose ) printf(" CALORIMETER - ERROR: RunInfo exited with non-zero status\n");
    code = sgnl;
    goto closeandexit;
  } else {
    sgnl = 0;
  };
  //
  // number of events in the file BEFORE the first event of our run
  //
  nobefrun = runinfo->GetFirstEntry();
  //
  // total number of events in the file 
  //
  totfileentries = runinfo->GetFileEntries();
  //
  // first file entry AFTER the last event of our run
  //
  noaftrun = runinfo->GetLastEntry() + 1;
  //
  // number of run to be processed
  //
  numbofrun = runinfo->GetNoRun();
  //
  // number of runs in the file
  //
  totnorun = runinfo->GetRunEntries();
  //
  // Does it contain already a Calorimeter branch? if so we are reprocessing data, if not we must create it.
  // 
  caloclone = (TTree*)file->Get("Calorimeter");
  //
  if ( !caloclone ){
    reproc = false;
    if ( run == 0 && verbose ) printf(" CALORIMETER - WARNING: you are reprocessing data but calorimeter tree does not exist!\n");
    if ( runinfo->IsReprocessing() && run != 0 && verbose ) printf(" CALORIMETER - WARNING: it seems you are not reprocessing data but calorimeter\n versioning information already exists in RunInfo.\n");
    //
  } else {
    //
    reproc = true;
    //
    if ( verbose ) printf("\n Preparing the pre-processing...\n");
    //
    if ( run == 0 ){
      //
      // if we are reprocessing everything we don't need to copy any old event and we can just create a new branch in the clone tree and jump steps 4/5/7.
      //
      if ( verbose ) printf("\n CALORIMETER - WARNING: Reprocessing all runs in the file\n");
      reprocall = true;
      //
    } else {
      //
      // we are reprocessing a single run
      //
      if ( verbose ) printf("\n CALORIMETER - WARNING: Reprocessing run number %u \n",run);
      reprocall = false;
      //
      //
      //
      tempfile = new TFile(tempname.str().c_str(),"RECREATE");
      tempcalo = caloclone->CloneTree(-1,"fast");
      tempcalo->SetName("Calorimeter-old");
      tempfile->Write();
      tempfile->Close();      
    };
    //
    // delete old tree
    //
    caloclone->Delete("all");
    //
    if ( verbose ) printf("\n ...done!\n");
    //
  };
  //
  // create calorimeter tree calo
  // 
  file->cd();
  calo = new TTree("Calorimeter-new","PAMELA Level2 calorimeter data");
  calo->Branch("CaloLevel2","CaloLevel2",&ca);
  //
  if ( reproc && !reprocall ){
    //
    // 
    //
    tempfile = new TFile(tempname.str().c_str(),"READ");
    caloclone = (TTree*)tempfile->Get("Calorimeter-old");
    caloclone->SetBranchAddress("CaloLevel2",&caclone);
    //      
    if ( nobefrun > 0 ){
      if ( verbose ) printf("\n Pre-processing: copying events from the old tree before the processed run\n");   
      if ( verbose ) printf(" Copying %u events in the file which are before the beginning of the run %u \n",nobefrun,run);
      if ( verbose ) printf(" Start copying at event number 0, end copying at event number %u \n",nobefrun);
      for (UInt_t j = 0; j < nobefrun; j++){
        //
        caloclone->GetEntry(j);   
        //
        // copy caclone to ca
        //
        memcpy(&ca,&caclone,sizeof(caclone));
        //
        // Fill entry in the new tree
        //
        calo->Fill();
        //
        ca->Clear();
        //
      };
      if ( verbose ) printf(" Finished successful copying!\n");
    };          
  };
  //
  // Get the list of run to be processed
  //
  runlist = runinfo->GetRunList();
  //
  // Loop over the run to be processed
  //
  for (UInt_t irun=0; irun < numbofrun; irun++){
    //
    badevent = 0;
    //
    idRun = runlist->At(irun);
    if ( verbose ) printf("\n\n\n ####################################################################### \n");
    if ( verbose ) printf("                    PROCESSING RUN NUMBER %u \n",idRun);
    if ( verbose ) printf(" ####################################################################### \n\n\n");
    //
    sgnl = runinfo->GetRunInfo(idRun);
    if ( sgnl ){
      if ( verbose ) printf(" CALORIMETER - ERROR: RunInfo exited with non-zero status\n");
      code = sgnl;
      goto closeandexit;
    } else {
      sgnl = 0;
    };
    id_reg_run = runinfo->ID_ROOT_L0;
    runheadtime = runinfo->RUNHEADER_TIME;
    runtrailtime = runinfo->RUNTRAILER_TIME;
    evfrom = runinfo->EV_FROM;
    evto = runinfo->EV_TO;
    //
    if ( id_reg_run == -1 ){
      if ( verbose ) printf("\n CALORIMETER - ERROR: no run with ID_RUN = %u \n\n Exiting... \n\n",idRun);
      code = -5;
      goto closeandexit;    
    };
    //
    // prepare the timesync for the db
    //
    if ( !dbc->IsConnected() ) throw -116;
    dbtime = new GL_TIMESYNC(runinfo->ID_ROOT_L0,"ID",dbc);
    //
    // Search in the DB the path and name of the LEVEL0 file to be processed.
    //
    if ( !dbc->IsConnected() ) throw -116;
    glroot->Query_GL_ROOT(runinfo->ID_ROOT_L0,dbc);
    //
    ftmpname.str("");
    ftmpname << glroot->PATH.Data() << "/";
    ftmpname << glroot->NAME.Data();
    fname = ftmpname.str().c_str();
    //
    // print out informations
    //
    totevent = runinfo->NEVENTS;
    if ( verbose ) printf("\n LEVEL0 data file: %s \n",fname.Data());
    if ( verbose ) printf(" RUN HEADER absolute time is:  %u \n",runheadtime);
    if ( verbose ) printf(" RUN TRAILER absolute time is: %u \n",runtrailtime);
    if ( verbose ) printf(" %i events to be processed for run %u: from %i to %i (reg entries)\n\n",totevent,idRun,evfrom,evfrom+totevent);
    //
    // Open Level0 file
    //
    l0File = new TFile(fname.Data());
    if ( !l0File ) {
      if ( verbose ) printf(" CALORIMETER - ERROR: problems opening  Level0 file\n");
      code = -6;
      goto closeandexit;
    };
    l0tr = (TTree*)l0File->Get("Physics");
    if ( !l0tr ) {
      if ( verbose ) printf(" CALORIMETER - ERROR: no Physics tree in Level0 file\n");
      l0File->Close();
      code = -7;
      goto closeandexit;
    };
    l0head = l0tr->GetBranch("Header");
    if ( !l0head ) {
      if ( verbose ) printf(" CALORIMETER - ERROR: no Header branch in Level0 tree\n");
      l0File->Close();
      code = -8;
      goto closeandexit;    
    };
    l0calo = l0tr->GetBranch("Calorimeter");    
    if ( !l0calo ) {
      if ( verbose ) printf(" CALORIMETER - ERROR: no Calorimeter branch in Level0 tree\n");
      l0File->Close();
      code = -103;
      goto closeandexit;
    };
    l0trig = l0tr->GetBranch("Trigger");
    if ( !l0trig ) {
      if ( verbose ) printf(" CALORIMETER - ERROR: no Trigger branch in Level0 tree\n");
      l0File->Close();
      code = -104;
      goto closeandexit;
    };
    //
    l0tr->SetBranchAddress("Trigger", &trig);
    l0tr->SetBranchAddress("Header", &eh);
    //
    // Construct the event object, look for the calibration which include the first header
    //  
    sgnl = 0;
    if ( verbose ) printf(" Check for calorimeter calibrations and initialize event object \n");
    if ( !dbc->IsConnected() ) throw -116;
    event->ProcessingInit(dbc,runheadtime,sgnl,l0tr,debug,verbose);
    if ( verbose ) printf("\n");
    if ( sgnl == 100 ) {
      code = sgnl;
      if ( verbose ) printf(" CALORIMETER - WARNING: run header not included in any calibration interval\n"); 
      sgnl = 0;
    };
    if ( sgnl ){
      l0File->Close();
      code = sgnl;
      goto closeandexit;
    };
    //
    qy.str("");
    //
    nevents = l0head->GetEntries();
    caloevents = l0calo->GetEntries();
    //
    if ( nevents < 1 ) {
      if ( verbose ) printf(" CALORIMETER - ERROR: Level0 file is empty\n\n");
      l0File->Close();
      code = -11;
      goto closeandexit;
    };
    //
    if ( evto > nevents-1 ) {
      if ( verbose ) printf(" CALORIMETER - ERROR: too few entries in the registry tree\n");
      l0File->Close();
      code = -12;
      goto closeandexit;
    };
    //
    // Check if we have to load parameter files
    //
    sgnl = 0;
    if ( !dbc->IsConnected() ) throw -116;
    sgnl = event->ChkParam(dbc,runheadtime); // calorimeter parameter files
    if ( sgnl < 0 ){
      code = sgnl;
      l0File->Close();
      goto closeandexit;
    };
    //
    if ( withtrk ){
      if ( trkpar1 || ( tttrkpar1 != 0 && tttrkpar1 < runheadtime ) ){
        trkpar1 = false;
        if ( !dbc->IsConnected() ) throw -116;
        Int_t glpar = q4->Query_GL_PARAM(runinfo->RUNHEADER_TIME,1,dbc);
        if ( glpar < 0 ){
          code = glpar;
          goto closeandexit;
        };
        tttrkpar1 = q4->TO_TIME;
        // ----------------------------
        // Read the magnetic field
        // ----------------------------
        if ( verbose ) printf(" Reading magnetic field maps at %s\n",(q4->PATH+q4->NAME).Data());
        trk->LoadField(q4->PATH+q4->NAME);
        if ( verbose ) printf("\n");
      };
    };
    //
    // run over all the events of the run
    //
    if ( verbose ) printf("\n Ready to start! \n\n Processed events: \n\n");
    //
    //
    for ( re = runinfo->EV_FROM; re < (runinfo->EV_FROM+runinfo->NEVENTS); re++){
      //
      if ( procev%1000 == 0 && procev > 0 && verbose ) printf(" %iK \n",procev/1000);   
      //
      l0head->GetEntry(re);
      //
      // absolute time of this event
      //
      ph = eh->GetPscuHeader();
      atime = dbtime->DBabsTime(ph->GetOrbitalTime());  
      //
      //
      //
      if ( re > caloevents-1 ){
        if ( verbose ) printf(" CALORIMETER - ERROR: no physics events with entry = %i in Level0 file\n",i);
        l0File->Close();
        code = -112;
        goto closeandexit;
      };
      //
      if ( atime > runtrailtime || atime < runheadtime  ) {
        if ( verbose ) printf(" CALORIMETER - WARNING: event at time outside the run time window, skipping it\n");
        jumped++;
        goto jumpev;
      };
      //
      // retrieve tracker informations
      //
      if ( !reprocall ){
        itr = nobefrun + (re - evfrom - jumped);
      } else {
        itr = runinfo->GetFirstEntry() + (re - evfrom - jumped);
      };
      //
      if ( withtrk ){
        if ( itr > nevtrkl2 ){  
          if ( verbose ) printf(" CALORIMETER - ERROR: no tracker events with entry = %i in Level2 file\n",itr);
          l0File->Close();
          code = -113;
          goto closeandexit; 
        };
        //
        trk->Clear();
        //
        tracker->GetEntry(itr); 
        //
      };
      //
      procev++;
      //
      // start processing
      //
      ca->Clear();
      //
      // determine who generate the trigger for this event (TOF, S4/PULSER, CALO) 
      //
      l0trig->GetEntry(re);
      S3 = 0;
      S2 = 0;
      S12 = 0;
      S11 = 0;
      S3 = trig->patterntrig[2];
      S2 = trig->patterntrig[3];
      S12 = trig->patterntrig[4];
      S11 = trig->patterntrig[5];
      if ( trig->patterntrig[0] ) tmptrigty = 2.;
      if ( S3 || S2 || S12 || S11 )  tmptrigty = 0.;    
      if ( trig->patterntrig[1] & (1<<0) || (!trig->patterntrig[0] && !S3 && !S2 && !S12 && !S11) ) tmptrigty = 1.;
      event->clevel2->trigty = tmptrigty;
      //
      // check if the calibration we are using is still good, if not load another calibration
      //
      sgnl = 0;
      sgnl = event->ChkCalib(dbc,atime);
      if ( sgnl < 0 ){
        code = sgnl;
        goto closeandexit;
      };
      if ( sgnl == 100 ){
        code = sgnl;
        if ( verbose ) printf(" CALORIMETER - WARNING: data not associated to any calibration interval\n"); 
        badevent++;
        sgnl = 0;      
      };
      //
      // do we have at least one track from the tracker? this check has been disabled
      //
      event->clevel1->good2 = 1;
      //
      //
      // Calibrate calorimeter event "re" and store output in the two structures that will be passed to fortran routine
      //
      event->Calibrate(re);      
      //
      // Here we have calibrated data, ready to be passed to the FORTRAN routine which will extract topological variables.
      //
      //
      // Calculate variables common to all tracks (qtot, nstrip, etc.)
      //
      event->GetCommonVar();
      //
      // Fill common variables
      //
      event->FillCommonVar(ca);
      //
      // Calculate variables related to tracks only if we have at least one track (from selftrigger and/or tracker)
      //
      ntrkentry = 0;
      //
      filled = false;
      //
      // Run over tracks (tracker or calorimeter )
      //
      if ( withtrk ){
        for(Int_t nt=0; nt < trk->ntrk(); nt++){  
          //
          event->clevel1->good2 = 1;
          //
          TrkTrack *ptt = trk->GetStoredTrack(nt); 
          //
          event->clevel1->trkchi2 = 0;
          //
          // Copy the alpha vector in the input structure
          //
          for (Int_t e = 0; e < 5 ; e++){
            event->clevel1->al_p[e][0] = ptt->al[e];
          };      
          //
          // Get tracker related variables for this track
          //
          event->GetTrkVar();
          //
          // Save tracker track sequence number
          //    
          event->trkseqno = nt;
          //
          // Copy values in the class ca from the structure clevel2
          //
          event->FillTrkVar(ca,ntrkentry);
          ntrkentry++;  
          filled = true;
          //
        }; // loop on all the tracks
      };
      // 
      // if no tracks found but there is the possibility to have a good track we should try to calculate anyway the track related variables using the calorimeter 
      // fit of the track (to be used for example when TRK is off due to any reason like IPM3/5 off).
      // here we make an event selection so it must be done very carefully...
      //
      // conditions are: 0) no track from the tracker 1) we have a track fit both in x and y 2) no problems with calo for this event 3) no selftrigger event
      //
      if ( trackanyway && !filled && event->clevel2->npcfit[0] >= 2 && event->clevel2->npcfit[1] >= 2 && event->clevel2->good != 0 && event->clevel2->trigty < 2. ){
        if ( debug ) printf(" Event with a track not fitted by the tracker at entry %i \n",itr);
        //
        // Disable "track mode" in the fortran routine
        //
        event->clevel1->good2 = 0;
        event->clevel1->riginput = rigdefault;
        if ( debug ) printf(" Using as default rigidity: %f \n",event->clevel1->riginput);
        //
        // We have a selftrigger event to analyze.
        //
        for (Int_t e = 0; e < 5 ; e++){
          event->clevel1->al_p[e][0] = 0.;
          event->clevel1->al_p[e][1] = 0.;
        };
        event->clevel1->trkchi2 = 0;
        //
        event->GetTrkVar();
        //
        // if we had no problem (clevel1->good2 = 0, NOTICE zero, not one in this mode!), fill and go on
        //
        if ( event->clevel1->good2 == 0 ) {
          //
          // In selftrigger mode the trkentry variable is set to -1
          //
          event->trkseqno = -3;
          //
          // Copy values in the class ca from the structure clevel2
          //
          event->FillTrkVar(ca,ntrkentry);
          ntrkentry++;
          filled = true;
          //
        } else {
          if ( verbose ) printf(" Selftrigger: problems with event at entry %i \n",itr);
        };
        //
      };
      //
      // Call high energy nuclei routine
      //
      if ( hZn && event->clevel2->trigty >= 2. ){
        if ( debug ) printf(" Calling selftrigger high energy nuclei routine for entry %i \n",itr);
        //
        // Disable "track mode" in the fortran routine
        //
        event->clevel1->good2 = 0;
        //
        // Set high energy nuclei flag to one
        // 
        event->clevel1->hzn = 1; 
        event->clevel1->riginput = rigdefault;
        //
        // We have a selftrigger event to analyze.
        //
        for (Int_t e = 0; e < 5 ; e++){
          event->clevel1->al_p[e][0] = 0.;
          event->clevel1->al_p[e][1] = 0.;
        };
        event->clevel1->trkchi2 = 0;
        //
        event->GetTrkVar();
        //
        // if we had no problem (clevel1->good2 = 0, NOTICE zero, not one in this mode!), fill and go on
        //
        if ( event->clevel1->good2 == 0 ) {
          //
          // In selftrigger mode the trkentry variable is set to -1
          //
          event->trkseqno = -2;
          //
          // Copy values in the class ca from the structure clevel2
          //
          event->FillTrkVar(ca,ntrkentry);
          ntrkentry++;
          filled = true;
          //
        } else {
          if ( verbose ) printf(" Selftrigger: problems with event at entry %i \n",itr);
        };
        //
      };
      //
      // self trigger event
      //
      if ( st && event->clevel2->trigty >= 2. ){
        if ( verbose ) printf(" Selftrigger event at entry %i \n",itr);
        //
        // Disable "track mode" in the fortran routine
        //
        event->clevel1->good2 = 0;
        //
        // disable high enery nuclei flag;
        //
        event->clevel1->hzn = 0; 
        //
        // We have a selftrigger event to analyze.
        //
        for (Int_t e = 0; e < 5 ; e++){
          event->clevel1->al_p[e][0] = 0.;
          event->clevel1->al_p[e][1] = 0.;
        };
        event->clevel1->trkchi2 = 0;
        //
        event->GetTrkVar();
        //
        // if we had no problem (clevel2->good = 0, NOTICE zero, not one in selftrigger mode!), fill and go on
        //
        if ( event->clevel1->good2 == 0 ) {
          //
          // In selftrigger mode the trkentry variable is set to -1
          //
          event->trkseqno = -1;
          //
          // Copy values in the class ca from the structure clevel2
          //
          event->FillTrkVar(ca,ntrkentry);
          ntrkentry++;
          filled = true;
          //
        } else {
          if ( verbose ) printf(" Selftrigger: problems with event at entry %i \n",itr);
        };
      };
      //
      if ( !filled ) badevent++;
      //
      // Clear structures used to communicate with fortran
      //
      event->ClearStructs();
      //
      // Fill the rootple 
      //
      calo->Fill();
      //    
      //      delete ca;
      //
    jumpev:
      debug = false;
      //
    };
    //
    if ( verbose ) printf("\n SUMMARY:\n");
    if ( verbose ) printf(" Total number of events: %i \n",totevent);
    if ( verbose ) printf(" Events with at least one track: %i \n",totevent-badevent);
    if ( verbose ) printf(" Events without tracks: %i \n",badevent);
    //
    if ( badevent == totevent ){
      if ( verbose ) printf("\n CALORIMETER - WARNING no tracks or good events in run %u \n",idRun);
      code = 101;    
    };
    //
    delete dbtime;
    //
    // Clear variables before processing another run  (needed to have always the same result when reprocessing data with the same software).
    //
    event->RunClose();
    if ( l0File ) l0File->Close();
    //
  }; // process all the runs
  //
  if ( verbose ) printf("\n Finished processing data \n");
  //
 closeandexit:
  //
  if ( !reprocall && reproc && code >= 0 ){
    if ( totfileentries > noaftrun ){
      if ( verbose ) printf("\n Post-processing: copying events from the old tree after the processed run\n");   
      if ( verbose ) printf(" Copying %u events in the file which are after the end of the run %u \n",(totfileentries-noaftrun),run);
      if ( verbose ) printf(" Start copying at event number %u end copying at event number %u \n",noaftrun,totfileentries);
      for (UInt_t j = noaftrun; j < totfileentries; j++ ){
        //
        // Get entry from old tree
        //
        caloclone->GetEntry(j);   
        //
        // copy caclone to ca
        //
        ca->Clear();
        //
        memcpy(&ca,&caclone,sizeof(caclone));
        //
        // Fill entry in the new tree
        //
        calo->Fill();
        //
      };
      if ( verbose ) printf(" Finished successful copying!\n");
    };
  };
  //
  // Case of no errors: close files, delete old tree(s), write and close level2 file 
  //
  if ( l0File ) l0File->Close();
  if ( tempfile ) tempfile->Close();            
  gSystem->Unlink(tempname.str().c_str());
  if ( tracker )  tracker->Delete();
  //
  if ( code < 0 ) printf("\n CALORIMETER - ERROR: an error occurred, try to save anyway...\n");
  if ( verbose ) printf("\n Writing and closing rootple\n");
  if ( runinfo ) runinfo->Close();    
  if ( calo ) calo->SetName("Calorimeter");    
  if ( file ){
    file->cd();
    file->Write();
  };
  if ( calo ) calo->Delete();
  //
  gSystem->Unlink(calofolder.str().c_str());
  //
  if ( ca ) delete ca;
  if ( caclone ) delete caclone;
  if ( trk ) delete trk;
  if ( q4 ) delete q4;
  if ( glroot ) delete glroot;
  if ( runinfo ) delete runinfo; 
  //
  // the end
  //
  if ( verbose ) printf("\n Exiting...\n");
  if ( code < 0 ) throw code;
  return(code);
}

