source: MEMPHYS/HEAD/applications/MEMPHYS_analysis.cxx @ 168

// AIDA analysis program over MEMPHYS.root file.

// AIDA :
#include <AIDA/AIDA.h>

// std :
#include <iostream>

//////////////////////////////////////////////////////////////////////////////
bool getXYZ(
 AIDA::ITuple& aParent
,int aColumn
,double& aX
,double& aY
,double& aZ
) 
//////////////////////////////////////////////////////////////////////////////
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
{
  AIDA::ITupleEntry* entry = (AIDA::ITupleEntry*)aParent.getObject(aColumn);
  if(!entry) return false;

  AIDA::ITuple* tuple = dynamic_cast<AIDA::ITuple*>(entry);
  if(!tuple) return false;

  tuple->start();
  if(!tuple->next()) return false;

  aX = tuple->getDouble(0);
  aY = tuple->getDouble(1);
  aZ = tuple->getDouble(2);

  return true;
}
//////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////
bool dump_tracks(
 AIDA::ITuple& aParent
) 
//////////////////////////////////////////////////////////////////////////////
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
{
  AIDA::ITupleEntry* entry = (AIDA::ITupleEntry*)aParent.getObject(8);
  if(!entry) return false;

  AIDA::ITuple* tuple = dynamic_cast<AIDA::ITuple*>(entry);
  if(!tuple) return false;

  //if(nPart!=tracks->rows()) {
  //  std::cout << "read: nPart / tracks mismatch " 
  //              << nPart << " " << tracks->rows() << std::endl;
  //  return false;
  //}

  tuple->start();
  int irow = 0;
  while(tuple->next()) {

    int pId = tuple->getInt(0);
    int parent = tuple->getInt(1);
    float timeStart = tuple->getFloat(2);

    double dx,dy,dz;
    if(!getXYZ(*tuple,3,dx,dy,dz)) return false;

    double mass = tuple->getDouble(4);
    double pTot = tuple->getDouble(5);
    double ETot = tuple->getDouble(6);

    double px,py,pz;
    if(!getXYZ(*tuple,7,px,py,pz)) return false;

    double start_x,start_y,start_z;
    if(!getXYZ(*tuple,8,start_x,start_y,start_z)) return false;

    double stop_x,stop_y,stop_z;
    if(!getXYZ(*tuple,9,stop_x,stop_y,stop_z)) return false;

    int startVol = tuple->getInt(10);
    int stopVol = tuple->getInt(11);

    std::cout << "----> Tk{"<<irow<<"}: " 
              << " pId " << pId
              << " parent " << parent
              << " creation time " << timeStart 
              << " Volumes " << startVol << " " << stopVol << "\n"
              << " Start Pos (" << start_x 
              << "," << start_y << "," << start_z << ")\n"
              << " Stop Pos (" << stop_x 
              << "," << stop_y << "," << stop_z << ")\n"
              << " dx,dy,dz " << dx << " " << dy << " " << dz << "\n"
              << " m " << mass
              << " ETot " << ETot
              << " pTot " << pTot
              << " px,py,pz " << px << " " << py << " " << pz << "\n"
              << std::endl;

    irow++;
  }
  return true;
}
//////////////////////////////////////////////////////////////////////////////
bool process_hits_times(
 AIDA::ITuple& aParent
,AIDA::IHistogram1D& aHisto
) 
//////////////////////////////////////////////////////////////////////////////
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
{
  AIDA::ITupleEntry* entry = (AIDA::ITupleEntry*)aParent.getObject(2);
  if(!entry) return false;

  AIDA::ITuple* tuple = dynamic_cast<AIDA::ITuple*>(entry);
  if(!tuple) return false;

  tuple->start();
  while(tuple->next()) {

    float time = tuple->getFloat(0);

    aHisto.fill(time);
  }
  return true;
}
//////////////////////////////////////////////////////////////////////////////
bool process_hits(
 AIDA::ITuple& aParent
,AIDA::IHistogram1D& aHisto
) 
//////////////////////////////////////////////////////////////////////////////
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
{
  AIDA::ITupleEntry* entry = (AIDA::ITupleEntry*)aParent.getObject(10);
  if(!entry) return false;

  AIDA::ITuple* tuple = dynamic_cast<AIDA::ITuple*>(entry);
  if(!tuple) return false;

  tuple->start();
  while(tuple->next()) {

    //int tubeId = tuple->getInt(0);
    //int totalPE = tuple->getInt(1);

    if(!process_hits_times(*tuple,aHisto)) return false;
  }
  return true;
}
//////////////////////////////////////////////////////////////////////////////
bool process_digits(
 AIDA::ITuple& aParent
,AIDA::IHistogram2D& aHisto
) 
//////////////////////////////////////////////////////////////////////////////
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
{
  AIDA::ITupleEntry* entry = (AIDA::ITupleEntry*)aParent.getObject(13);
  if(!entry) return false;

  AIDA::ITuple* tuple = dynamic_cast<AIDA::ITuple*>(entry);
  if(!tuple) return false;

  tuple->start();
  while(tuple->next()) {

    //int tubeId = tuple->getInt(0);
    double pe = tuple->getDouble(1);
    double time = tuple->getDouble(2);
    //printf("debug : ++++ : %g %g\n",time,pe);

    aHisto.fill(time,pe);
  }
  return true;
}
//////////////////////////////////////////////////////////////////////////////
bool plot(
 AIDA::IAnalysisFactory& aAIDA
,AIDA::IHistogram1D& aHisto1D
,AIDA::IHistogram2D& aHisto2D
) 
//////////////////////////////////////////////////////////////////////////////
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
{
  std::string header = "MEMPHYS_analysis : ";

  AIDA::IPlotterFactory* plotterFactory = 
    aAIDA.createPlotterFactory(0,0);
  if(!plotterFactory) {
    std::cout << header << "can't get a PlotterFactory." << std::endl;
    return false;
  }
    
  AIDA::IPlotter* plotter = plotterFactory->create();
  if(!plotter) {
    std::cout << header << "can't get a plotter." << std::endl;
    return false;
  }
    
  plotter->clearRegions();
  plotter->createRegions(1,2,0);

  plotter->setCurrentRegionNumber(0);
  plotter->currentRegion().plot(aHisto1D);

  plotter->setCurrentRegionNumber(1);
  plotter->currentRegion().plot(aHisto2D);


  plotter->show();
  plotter->interact();

  delete plotter;
  delete plotterFactory;

  return true;
}
//////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////
int main(
 int //aArgc
,char** //aArgv
) 
//////////////////////////////////////////////////////////////////////////////
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
{
  //std::string arg1 = aArgc>=2?aArgv[1]:"";

  std::string program = "MEMPHYS_analysis";
  std::string header = program + " : ";

  AIDA::IAnalysisFactory* aida = AIDA_createAnalysisFactory();
  if(!aida) {
    std::cout << header << "AIDA not found." << std::endl;
    return 1;
  }

  AIDA::ITreeFactory* treeFactory = aida->createTreeFactory();
  if(!treeFactory) {
    std::cout << header << "can't get a TreeFactory." << std::endl;
    return 1;
  }

  ////////////////////////////////////////////////////////
  // Create histograms in memory tree ////////////////////
  ////////////////////////////////////////////////////////
  AIDA::ITree* memory = treeFactory->create();
  if(!memory) {
    std::cout << header << "can't get memory tree." << std::endl;
    return 1;
  }
  AIDA::IHistogramFactory* histogramFactory = 
    aida->createHistogramFactory(*memory);
  if(!histogramFactory) {
    std::cout << header << "can't get an histogram factory." << std::endl;
    return 1;
  }
  
  AIDA::IHistogram1D* hits_times = 
    histogramFactory->createHistogram1D("hits_times","Hits times",100,0,3000);
  if(!hits_times) {
    std::cout << header 
              << "can't create histogram : time." 
              << std::endl;
    return 1;
  }

  AIDA::IHistogram2D* digits_time_pe = 
    histogramFactory->createHistogram2D("digits_pe_time","Digits PE time",
                                        100,0,3000,100,0,10);
  if(!digits_time_pe) {
    std::cout << header 
              << "can't create histogram : digits_time_pe." 
              << std::endl;
    return 1;
  }

  delete histogramFactory;

  ////////////////////////////////////////////////////////
  /// Read data //////////////////////////////////////////
  ////////////////////////////////////////////////////////

  AIDA::ITree* tree = treeFactory->create("MEMPHYS.root","root",true,false);
  if(!tree) {
    std::cout << header << "can't open data file." << std::endl;
    return 1;
  }

  AIDA::IManagedObject* object = tree->find("Event");
  if(!object) {
    std::cout << header 
              << "object Event not found in tree." 
              << std::endl;
    return 1;
  }
  AIDA::ITuple* tuple = dynamic_cast<AIDA::ITuple*>(object);
  if(!tuple) {
    std::cout << header << "object not an AIDA::ITuple." << std::endl;
    return 1;
  }

  int coln = tuple->columns();
  for(int index=0;index<coln;index++) {
    std::cout << header 
              << " icol = " << index
              << ", label = " << tuple->columnName(index) 
              << ", type = " << tuple->columnType(index) 
              << std::endl;
  }

  std::cout << header << "rows = " << tuple->rows() << std::endl;

  int nentries = tuple->rows();
  //int nentries = 100000;
  //int nentries = 40;
  std::cout << header 
            << "traitements de " << nentries << " entrees" 
            << std::endl;

  int irow = 0;
  tuple->start();
  while(tuple->next() && (irow<nentries)) {

    /*
    int eventId = tuple->getInt(0);
    //int inputEvtId = tuple->getInt(1);
    //int interMode = tuple->getInt(2);
    //int vtxVol = tuple->getInt(3);

    int nPart = tuple->getInt(5);
    //int leptonIndex = tuple->getInt(6);
    //int protonIndex = tuple->getInt(7);

    int nHits = tuple->getInt(9);
    int nDigits = tuple->getInt(11);
    double sumPE = tuple->getDouble(12);

    std::cout << ">>>>>>>>>>>>> Event{" << irow << "}: "
              << " evt Id " << eventId 
    //        << " evt Input Id " << inputEvtId
    //        << "\n interaction mode " << interMode
    //        << " start in volume " << vtxVol << "\n"
              <<" #tracks: " << nPart
              <<" #hits: " << nHits
              <<" #digits: " << nDigits
              <<" sumPE " << sumPE
              << std::endl;
    */

    //if(!dump_tracks(*tuple)) return 1;

    if(!process_hits(*tuple,*hits_times)) return 1;
    if(!process_digits(*tuple,*digits_time_pe)) return 1;

    irow++;
  }

  delete treeFactory;
  delete tree;

  plot(*aida,*hits_times,*digits_time_pe);

  delete aida;

  return 0;
}