/data3/calcul/jacquem/working_dir/Micromegas/micromegasFrameWork/src/reconstruction.cpp

Go to the documentation of this file.

/* @version $Revision: 1945 $ * @modifiedby $Author: jacquem $ * @lastmodified $Date: 2012-11-20 10:37:04 +0100 (Tue, 20 Nov 2012) $ */

#include "geometry/Chamber.hh"
#include "geometry/Board.hh"
#include "geometry/Chip.hh"
#include "geometry/Channel.hh"
#include "geometry/Dif.hh"
#include "geometry/Detector.hh"

#include "event/Run.hh"
#include "event/Event.hh"

#include "parser/AcquisitionParser.hh"
#include "parser/Centaure.hh"
#include "parser/Hardroc1Reader.hh"
#include "parser/Hardroc2XdaqReader.hh"
#include "parser/Hardroc2LabviewReader.hh"
#include "parser/MicrorocOldLabviewReader.hh"
#include "parser/MicrorocLabviewReader.hh"
#include "parser/MicrorocXDaqReader.hh"
#include "parser/DifSynchroReader.hh"
#include "parser/DiracReader.hh"
#include "parser/DiracLabview.hh"
#include "parser/CalibHR1Parser.hh"
#include "parser/CalibHR2Parser.hh"
#include "parser/CalibMicrorocParser.hh"
#include "parser/TestMicrorocParser.hh"

#include "xml/XMLTool.hh"

#include "tools/Log.hh"
#include "tools/MicroException.hh"
#include "tools/Arguments.hh"

#include "MTEvent.hh"
#include "MTRun.hh"
#include "slowControl/SlowControlManager.hh"
#include "slowControl/SlowControlEntry.hh"

#include "mysql/Mysql.hh"

#include <iostream>
#include <iomanip>
#include <sstream>
#include <iostream>

#include <stdlib.h>


#include "TFile.h"
#include "TTree.h"



#include "IO/LCWriter.h"
#include "IMPL/LCRunHeaderImpl.h"
#include "IMPL/LCEventImpl.h"
#include "IMPL/LCTOOLS.h"

#include "lcio/MLcio.hh"

using namespace lcio;
using namespace std;

int main(int argc, char **argv)
{
  string steerName;


  // set verbose level
  FILELog::ReportingLevel() = FILELog::FromString(INFO);

  ui32 slowcontrolNotFound = 0;
  ui32 eventNotValid = 0;

//----- Control usage and option
  Arguments arg(argc,argv,"reconstruction steerfile.xml type(root|lcio) [-r runId]  [-l nbEvent] -dif[id1,id2,...] [-sql 1]  [-forceanalog 1] [-q 1] [-batch 1] [-skip nbEvent]  [-v(ERROR || WARNING ||  INFO  || DEBUG || DEBUG1)]");
  
  if ( arg.getNbOfArgs() != 2) { //  steerfile.xml type
    arg.usage();
    exit(1);
  }
  // Get command parameter : here steering file
  steerName.assign(arg.getArgument(1));
  
  string outputType = arg.getArgument(2);

  bool outputLcio = true;
  if ( outputType.compare("root") == 0 )
  {   
    outputLcio = false;
  }
  
  
  bool batch = false;
  if ( arg.getOption("-batch").size() != 0 )
  {
    if ( atoi(arg.getOption("-batch").c_str()) == 1) 
    {
      cout << "Batch mode" << endl;
      batch = true;
    }
  }

  if ( arg.getOption("-v").size() != 0 )
  {
    cout << "VERBOSE LEVEL: " << arg.getOption("-v") << endl;
    FILELog::ReportingLevel() = FILELog::FromString(arg.getOption("-v"));
  }

  bool sql = false;
  if ( arg.getOption("-sql").size() != 0 )
  {
    cout << "Mysql database will be updated" << endl;
    sql = true;
  }

  ui32 limitEvent = 0;
  string sLimitEvent;
  bool limit = false;
  if ( arg.getOption("-l").size() != 0 ) 
  {
    limitEvent = atoi(arg.getOption("-l").c_str());
    sLimitEvent = arg.getOption("-l");
    limit = true;
    FILE_LOG(logINFO)  <<"Number of reconstruct event is limited to " << limitEvent << endl;
  } 

  ui32 skipEvent = 0;
  //string sSkipEvent;
  bool skip = false;
  if ( arg.getOption("-skip").size() != 0 ) 
  {
    skipEvent = atoi(arg.getOption("-skip").c_str());
    skip = true;
    FILE_LOG(logINFO)  <<"Skip " << skipEvent << " events." << endl;
  } 

  set<unsigned int>*difList = NULL;
  if ( arg.getOption("-dif").size() != 0 )
  {
    difList = new set<unsigned int>();
    string sdifList = arg.getOption("-dif");
    string difId;
    stringstream stream(sdifList);
    cout  <<"Create hit only for following dif: " ;
    while( getline(stream, difId, ',') )
    {
      difList->insert(atoi(difId.c_str()));
      cout << atoi(difId.c_str()) << ", ";
    }
    cout << endl;
  }

  // Add by Renaud to limit warnings during microroc qualifaction (no database is used)
  bool is_qualification = false;
  if ( arg.getOption("-q").size() != 0 ) 
  {
        is_qualification = true;
  }
  
  /* alow to create analog hits event if no digital hit correspond. Warning hit will have no bcid information
                             // http://lappwiki01.in2p3.fr/Linear_Collider_LAPP/index.php/MicromegasFramework#Reconstruction_de_la_lecture_analogique_avec_les_Microroc
  */
  bool forceAnalog = false;
  if ( arg.getOption("-forceanalog").size() != 0 ) 
  {
        forceAnalog = true;
  }

  string svnVersion = "0";
  // test if steer file exist
  FILE *file = fopen(steerName.c_str(), "r");
  if(file == NULL){
    FILE_LOG(logERROR) << "Steer file ["<< steerName.c_str() << "] does not exist" << endl;
    exit(1);
    fclose(file);
  }
  else
  {
    // create Run, set detector and steerDesc from steer XML file
    SteerDesc steerDesc;
    XMLTool xml;
    xml.parse(steerDesc, steerName); // parse steer file and fill steerDesc

    fclose(file);

    string slowControlPath = steerDesc.getSlowControlPath();
    SlowControlManager sc;
    sc.build(steerDesc);
    sc.readDB(slowControlPath);

    FILE_LOG(logINFO) <<  "slowControl entries:[" << sc.getNbEntries() << "]" << endl;

    TFile *hfile = NULL;
    MTEvent *rootEvt = NULL;

    LCWriter* lcWrt = NULL;
    LCRunHeaderImpl* runHdr = NULL ;

    if ( outputLcio )
    {
      runHdr = new LCRunHeaderImpl() ;

      string detName("micromegas_m2")  ;
      runHdr->setDetectorName( detName ) ;

      string ecalName("m2") ;
      runHdr->addActiveSubdetector( ecalName ) ;

      string telName("telecsope") ;
      runHdr->addActiveSubdetector( telName ) ;
      // add some parameters to the run header
      runHdr->parameters().setValue( "cerate lcio for micromegas" , "lcioTest.cc" ) ;
      IntVec iv(3) ;
      iv[0] = 1 ;
      iv[1] = 2 ;
      iv[2] = 3 ;
      runHdr->parameters().setValues( "SomeIndices" , iv ) ;
      runHdr->parameters().setValues( "SomeOtherIndices" , iv ) ;
    }
    else 
    {
      string rootFile = "default.root";
      if ( steerDesc.outputFile.length() != 0)  { rootFile = steerDesc.outputFile; }
      hfile = new TFile(rootFile.data(),"RECREATE","Micromegas LAPP");

      // Start ROOT tree filling
      rootEvt = new MTEvent();   // MTEvent is a Root base class which can be write to a ROOT File.
    }
    unsigned int totalNbEvents = 0;
    unsigned int totalNbHits = 0;
    unsigned int totalNbEventPerRawFile = 0;
    UInt_t lastEventId = 0;
    // reconstruct event from input data files
    FILE_LOG(logINFO) << "Start reconstruction for " << steerName  << endl;
    FILE_LOG(logINFO) << "Output file " << steerDesc.outputFile  <<  endl;
    FILE_LOG(logINFO) <<  "Output format " << outputType << "." << endl;

    Detector detector;
    if ( ! detector.build(steerDesc) == steerDesc.chambers.size() )
    {
      FILE_LOG(logERROR) << "Geometry Error. Please check your XML file."  << endl;
      exit(-1);
    }
    Run run(detector);
    svnVersion = run.getSvnrev();
    if ( ! steerDesc.setRun(run)) // Set detector with xml file content
    { exit(1); }

    // create the proper parser depending on file type
    const std::vector<input_info_t>* inputFiles = &(steerDesc.inputFiles);

    ui16 nbFile =  inputFiles->size();
    for (int fileIndex = 0; fileIndex < nbFile; fileIndex++)
    {
      const input_info_t &inputInfo = inputFiles->at(fileIndex);
      std::string currentType = inputInfo.type;
      FILE *inputFile = fopen(inputInfo.path.data(),"rb");
      if ( inputFile == NULL)
      {
        FILE_LOG(logERROR) << "Could not open " << inputInfo.path << endl;
        exit(1);
        continue;
      }


      FILE_LOG(logINFO) << "Open acquisition file: " << inputInfo.path << " of type:" << currentType << endl;
      AcquisitionParser *parser;

      if (currentType == Centaure::type())
      {
        parser = new Centaure(run, inputFile, lastEventId);
      }
      else if (currentType == Hardroc1Reader::type())
      {
        parser = new Hardroc1Reader(run, inputFile, lastEventId);
      }
      else if (currentType == Hardroc2XdaqReader::type())
      {
        parser = new Hardroc2XdaqReader(run, inputFile, lastEventId);
      }
      else if (currentType == Hardroc2LabviewReader::type())
      {
        parser = new Hardroc2LabviewReader(run, inputFile, lastEventId);
      }
      else if (currentType == CalibHR1Parser::type())
      {
        parser = new CalibHR1Parser(run, inputFile, lastEventId);
      }
      else if (currentType == CalibHR2Parser::type())
      {
        parser = new CalibHR2Parser(run, inputFile, lastEventId);
      }
      else if (currentType == MicrorocOldLabviewReader::type())
      {
        parser = new MicrorocOldLabviewReader(run, inputFile, lastEventId);
      }
      else if (currentType == CalibMicrorocParser::type())
      {
        parser = new CalibMicrorocParser(run, inputFile, lastEventId);
      }
      else if (currentType == TestMicrorocParser::type())
      {
        parser = new TestMicrorocParser(run, inputFile, lastEventId);
      }
      else if (currentType == DiracReader::type())
      {
        parser = new DiracReader(run, inputFile, lastEventId);
      }
      else if (currentType == DiracLabview::type())
      {
        parser = new DiracLabview(run, inputFile, lastEventId);
      }
      else if (currentType == DifSynchroReader::type())
      {
        parser = new DifSynchroReader(run, inputFile, lastEventId);
      }
      else if (currentType == MicrorocLabviewReader::type())
      {
        parser = new MicrorocLabviewReader(run, inputFile, lastEventId,forceAnalog, difList);
        if ( skip)  { parser->setSkipEvent(skipEvent) ; } 
      }
      else if (currentType == MicrorocXDaqReader::type())
      {
        parser = new MicrorocXDaqReader(run, inputInfo.path.data(), lastEventId,forceAnalog, difList);
        if ( skip)  { parser->setSkipEvent(skipEvent) ; } 
      }
      else if (currentType == "xdaqHR2")
      {
        parser = new MicrorocXDaqReader(run, inputInfo.path.data(), lastEventId,difList);
      }

      else
      { // type non pris en charge
      FILE_LOG(logERROR) << "NOT SUPPORTED ---- type:" << currentType << endl;
      fclose(inputFile);
      continue;
      }

      totalNbEventPerRawFile = 0;

      unsigned int configNum = 1;
      Event *event = new Event(run);

      ui16 parserResult = parser->getNextEvent(*event);
      while (parserResult == NEW_CONFIG)
      {  
        size_t  found=inputInfo.path.find_last_of("/\\");
        string name = inputInfo.path.substr(found+1);
        string s = "chaine";
        stringstream buffer;
        buffer << name << "_" << configNum ;


        TTree mt(buffer.str().c_str(),"MicroMegas event");

        if (  outputLcio )
        {
        // create sio writer
        cout << " runHdr->setRunNumber[" << inputInfo.runId << "]" << endl;
        runHdr->setRunNumber(inputInfo.runId ) ;
        lcWrt = LCFactory::getInstance()->createLCWriter()  ;
        string lcioFile = buffer.str();
        lcWrt->open( lcioFile , LCIO::WRITE_NEW )  ;
        }
        else
        {
        //              FILE_LOG(logINFO) << "Create new TTRee:["<< buffer.str() << "]" << endl;
        mt.SetMaxTreeSize(40000000000); // 40 GB
        mt.Branch("MTEvent",&rootEvt,16000,2);// second argument is class of evt!!!
        }

        // now parse the whole file
        // now fill the run structure
        // mtRun = run;
        configNum++;

        parserResult = parser->getNextEvent(*event);

        run.setRawDataFilename(inputInfo.path);
//        MTRun mtRun(run);

        while (parserResult == EVENT_CORRECT || parserResult == EVENT_ERROR  ) // && ( !limit || totalNbEventPerRawFile < limitEvent  ))
        {
          if ( parserResult == EVENT_CORRECT )
          {
            try
            {
              if ( sc.getNbEntries() > 0 )
              {
              const double evTimeStamp = event->getTimeStamp();
                try
                {
                  const SlowControlEntry* slowcontrolEntry = sc.getEntry(evTimeStamp);
                  event->setSlowControlParam(slowcontrolEntry);
                }
                catch ( MicroException e)
                {
                  FILE_LOG(logWARNING) << e.getMessage() << endl;
                  slowcontrolNotFound++;
                }
              }
              if ( ! outputLcio )
              {
                *rootEvt=*event;   // ! surcharge operator=
              }
              totalNbHits =  totalNbHits + event->getChannelHitSize();

              if ( event->getId() % 1== 0)
              {
                if ( batch )
                {
                  FILE_LOG(logINFO  ) << "Write Event[" << event->getId() << "] total number of reconstruted event' hits[" << event->getChannelHitSize() << "] valid[" << event->getValidFlag() << "]" << endl;
                }
                else
                {  
                  if ( ! event->isValid()  ) 
                  {
                    FILE_LOG(logINFO  ) << "Write Event[" << event->getId()-1 << "] total number of reconstruted event' hits[" << event->getChannelHitSize() << "] , FLAG[" << hex  << event->getValidFlag() <<  "] ->["  << event->getValidFlagString() << "]" << dec <<endl;
                  }
                  else
                  {
                    FILE_LOG(logINFO  ) << "Write Event[" << event->getId()-1 << "] total number of reconstruted event' hits[" << event->getChannelHitSize()  << "\t\t\t\t   \r" << flush;
                  }
                }
                //if ( event->getValidFlag() != VALID ) { eventNotValid++; }
                if ( ! event->isValid()) { eventNotValid++; }
                
              }

              lastEventId = event->getId();

              if (  outputLcio )
              {
                // write the event to the file
                LCEventImpl* lcioEvtImpl =  new LCEventImpl();
                MLcio::setLCEventImp(*lcioEvtImpl,*event);
                lcWrt->writeEvent( lcioEvtImpl ) ;
                // dump the event to the screen
                //LCTOOLS::dumpEvent( lcioEvtImpl ) ;
                // ------------ IMPORTANT ------------- !
                // we created the event so we need to delete it ...
                delete lcioEvtImpl ;
              }

              else 
              {
                mt.Fill(); // Fill the TTree
                rootEvt->Clear(); // free memory
              }
              totalNbEvents++;
              totalNbEventPerRawFile++;

            }
            catch (MicroException e) { FILE_LOG(logWARNING) << e.getMessage() << endl; }
            // Save all objects in this file
            delete event;
          } // END  if ( parserResult == EVENT_CORRECT)
          if ( !limit || totalNbEventPerRawFile < limitEvent  )
          {
            event = new Event(run);
            parserResult = parser->getNextEvent(*event);
          }
          else
          {
            parserResult = LAST_FILE_EVENT;
          }
        } // end for Events
        MTRun mtRun(run);

        if (  outputLcio )
        {
          lcWrt->close() ;
          delete lcWrt;
        }
        else if ( sql )
        {
          // get Run fronm database
          Mysql mysql;
          const input_info_t &inputInfo = inputFiles->at(0);
          string full = inputInfo.path.data();
          size_t last_slash = full.find_last_of("/")  +1 ;
          string name = full.substr(last_slash);
          string path = full.substr(0,last_slash-1);
          ui32 runId = mysql.getRunIdFromRawFile(path,name);
          if ( runId != UNDEFINED_ID )
          { 
            mtRun.SetRunId(runId);
            FILE_LOG(logINFO) << "Set Run Id from database[" << runId << "]"  << endl;
          }
        }
        mt.GetUserInfo()->Add(&mtRun);  // store run information in a TLIst a mt TTree
        TFile *curFile =mt.GetCurrentFile() ;
        curFile->Write("", TObject::kOverwrite);

        run.setDataFormat(0);
      }

      // free the memory and close the file
      delete parser;
      fclose(inputFile);
      FILE_LOG(logINFO) << endl << "-- " << inputInfo.path.data() <<
                        " Number of reconstructed Events:[" << totalNbEventPerRawFile << "]" << endl <<  endl;
    } // end for file index



    FILE_LOG(logINFO) << endl << "----- Total number of reconstructed Events:[" << totalNbEvents << "]" << endl;
    FILE_LOG(logINFO) << endl << "----- Total number of Hits:[" << totalNbHits << "]" << endl;

    if ( sql)
    {
      Mysql mysql;
      // Add reconstruction software in  DB if not already exist
      FILE_LOG(logINFO) << " mysql.add reconstruction software to DB (" << svnVersion << "]"  << endl;
      int softId = mysql.addReconstructionSoftware(svnVersion);

      // Add reconstruct file
      string fullName = steerDesc.outputFile ;
      size_t last_slash = fullName.find_last_of("/")  +1 ;
      string name = fullName.substr(last_slash);
      string path = fullName.substr(0,last_slash);
      string comment = "No comment.";
      if ( limitEvent != 0 ) { comment = "Reconstrcuted event limited to " + sLimitEvent + "." ; }
      i16 rebuildId = mysql.addRebuildFile(softId,path, name, true, totalNbEvents , totalNbHits, true, comment );

      // Get Raw file database Id if they are already existing
      for (int fileIndex = 0; fileIndex < inputFiles->size(); fileIndex++)
      {
        const input_info_t &inputInfo = inputFiles->at(fileIndex);
        string full = inputInfo.path.data();
        size_t last_slash = full.find_last_of("/")  +1 ;
        string name = full.substr(last_slash);
        string path = full.substr(0,last_slash-1);
        try 
        { 
          ui32 rawId = mysql.getRawFileId(path,name);
          mysql.connectRebuildAndRawFile(rawId,rebuildId);
        }
        catch ( MicroException e ) 
        {
          if (is_qualification == false)
          FILE_LOG(logWARNING) << "Raw file[" << path + name << "] not registerd in database." << endl;
        }
      }
    }
    /*
    The default value for AutoSave is to save a copy of the Tree header
    as soon as you have written 100 Mbytes (uncompressed) to the file.
    When AutoSave is called, a new key is created with the current header.
    If the creation of the key is successful, the previous key is deleted.
    At the end of the job, when you do tree.Write, a new key is
    always created, keeping the previous cycle of the same key, if any.
    You can use the option kOverwrite when calling the final tree.Write
    if you want to keep one single key.
    */
    if ( ! outputLcio )
    {
      hfile->Close();
    }
    FILE_LOG(logINFO) << "Number of slowcontrol data not found:" << slowcontrolNotFound << endl;
    FILE_LOG(logINFO) << "Number of non valid event: " << eventNotValid  << endl;
  }
  
  return 0;
}

---