// AIDA analysis program over ELYSE.root file. // AIDA : #include #include // std : #include // Lib : #include ////////////////////////////////////////////////////////////////////////////// AIDA::ITuple* cast_Tuple( AIDA::ITupleEntry* aEntry ) ////////////////////////////////////////////////////////////////////////////// //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!// { //FIXME : CINT can't handle dynamic_cast : //return (AIDA::ITuple*)aEntry->cast("AIDA::ITuple"); return dynamic_cast(aEntry); } ////////////////////////////////////////////////////////////////////////////// bool get_XYZ( 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 = cast_Tuple(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 = cast_Tuple(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(!get_XYZ(*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(!get_XYZ(*tuple,7,px,py,pz)) return false; double start_x,start_y,start_z; if(!get_XYZ(*tuple,8,start_x,start_y,start_z)) return false; double stop_x,stop_y,stop_z; if(!get_XYZ(*tuple,9,stop_x,stop_y,stop_z)) return false; int startVol = tuple->getInt(10); int stopVol = tuple->getInt(11); std::cout << "----> Tk{"<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 = cast_Tuple(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 = cast_Tuple(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; } ////////////////////////////////////////////////////////////////////////////// void set_region_style( AIDA::IPlotterRegion& aRegion ) ////////////////////////////////////////////////////////////////////////////// // Taken from G4SPL_analysis.cxx. //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!// { //FIXME : have to do the below with AIDA styles. // ROOT is in NDC, then we take a plotter with width = height = 1 // PAW : //float XSIZ = 20.0F; //float XLAB = 1.4F; //float YSIZ = 20.0F; //float YLAB = 1.4F; float XLAB = 1.4F/20.0F; //0.07 //x distance of y title to data frame. float YLAB = 0.8F/20.0F; //0.04 //y distance of x title to data frame. YLAB = 0.05F; //FIXME : cooking. // To have a good matching with ROOT for text size : //double majic = 0.014/0.027; double majic = 0.6; std::string s = Lib::smanip::tostring(majic); //FIXME : have methods for the below : aRegion.style().setParameter("textScale",s); aRegion.style().xAxisStyle().setParameter("magStyle.scale",s); aRegion.style().yAxisStyle().setParameter("magStyle.scale",s); aRegion.style().xAxisStyle().tickLabelStyle().setParameter("scale",s); aRegion.style().xAxisStyle().labelStyle().setParameter("scale",s); aRegion.style().yAxisStyle().tickLabelStyle().setParameter("scale",s); aRegion.style().yAxisStyle().labelStyle().setParameter("scale",s); // ROOT def margins 0.1. SoPlotter def 0.1. aRegion.layout().setParameter("rightMargin",0.1); aRegion.layout().setParameter("topMargin",0.1); //gStyle->SetPadBottomMargin(0.15); aRegion.layout().setParameter("leftMargin",0.15); //gStyle->SetPadLeftMargin(0.15); aRegion.layout().setParameter("bottomMargin",0.15); aRegion.style().setParameter("superposeBins","TRUE"); aRegion.setParameter("plotter.wallStyle.visible","FALSE"); aRegion.setParameter("plotter.gridStyle.visible","FALSE"); //std::string font = "Hershey"; //PAW font. //std::string font = "TTF/arialbd"; //Default ROOT font 62. //std::string font = "TTF/couri"; std::string font = "TTF/times"; //= ROOT 132 ? //std::string smoothing = "FALSE"; std::string smoothing = "TRUE"; // X axis : aRegion.style().xAxisStyle().tickLabelStyle().setFont(font); aRegion.style().xAxisStyle().tickLabelStyle().setParameter("smoothing",smoothing); aRegion.style().xAxisStyle().labelStyle().setFont(font); aRegion.style().xAxisStyle().labelStyle().setParameter("smoothing",smoothing); aRegion.setParameter("plotter.xAxis.magStyle.fontName",font); aRegion.setParameter("plotter.xAxis.magStyle.smoothing",smoothing); aRegion.style().xAxisStyle().lineStyle().setThickness(2); aRegion.setParameter("plotter.xAxis.ticksStyle.width","2"); aRegion.style().setParameter("topAxisVisible","TRUE"); // Set hplot tick modeling : aRegion.style().xAxisStyle().setParameter("modeling","hplot"); aRegion.style().xAxisStyle().setParameter("divisions","505"); // Y axis : //aRegion.style().setParameter("yAxisLogScale","TRUE"); aRegion.style().yAxisStyle().tickLabelStyle().setFont(font); aRegion.style().yAxisStyle().tickLabelStyle().setParameter("smoothing",smoothing); aRegion.style().yAxisStyle().labelStyle().setFont(font); aRegion.style().yAxisStyle().labelStyle().setParameter("smoothing",smoothing); aRegion.setParameter("plotter.yAxis.magStyle.fontName",font); aRegion.setParameter("plotter.yAxis.magStyle.smoothing",smoothing); aRegion.style().yAxisStyle().lineStyle().setThickness(2); aRegion.setParameter("plotter.yAxis.ticksStyle.width","2"); //gStyle->SetTitleSize(0.06,"XYZ"); //ROOT def 0.04. SoPlotter : 0.014 aRegion.setParameter("plotter.xAxis.titleHeight","0.06"); aRegion.setParameter("plotter.yAxis.titleHeight","0.06"); //gStyle->SetLabelOffset(0.01,"Y"); //ROOT def 0.005. SoPlotter def 0.02 aRegion.setParameter("plotter.yAxis.labelToAxis","0.01"); //gStyle->SetTitleOffset(1.1,"Y"); aRegion.setParameter("plotter.yAxis.titleToAxis", Lib::smanip::tostring(1.1*XLAB)); aRegion.setParameter("plotter.xAxis.titleToAxis", Lib::smanip::tostring(1.1*YLAB)); //gStyle->SetLabelSize(0.05,"XYZ"); //ROOT def 0.04. SoPlotter def 0.014. aRegion.setParameter("plotter.yAxis.labelHeight","0.05"); aRegion.setParameter("plotter.xAxis.labelHeight","0.05"); aRegion.style().setParameter("rightAxisVisible","TRUE"); // Set hplot tick modeling : aRegion.style().yAxisStyle().setParameter("modeling","hplot"); aRegion.style().yAxisStyle().setParameter("divisions","505"); // title : aRegion.style().setParameter("titleHeight","0.04"); aRegion.style().setParameter("titleToAxis","0.02"); aRegion.style().titleStyle().textStyle().setFont(font); aRegion.style().titleStyle().textStyle().setParameter("smoothing",smoothing); aRegion.style().titleStyle().textStyle().setColor("0 0 1"); // legend box : aRegion.setParameter("legendRegionVisible","TRUE"); aRegion.setParameter("legendRegion.viewportRegion.backgroundColor","1 1 0"); aRegion.setParameter("legendRegionSize","0.15 0.16"); aRegion.setParameter("legendRegionOrigin","0.1 0.1"); aRegion.setParameter("legendRegion.horizontalMargin","2"); aRegion.setParameter("legendRegion.verticalMargin","2"); // Frame : aRegion.setParameter("plotter.innerFrameStyle.lineWidth","2"); } ////////////////////////////////////////////////////////////////////////////// bool plot( AIDA::IAnalysisFactory& aAIDA ,AIDA::IHistogram1D& aHisto1D ,AIDA::IHistogram2D& aHisto2D ) ////////////////////////////////////////////////////////////////////////////// //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!// { std::string header = "ELYSE_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); AIDA::IPlotterRegion& region = plotter->currentRegion(); set_region_style(region); region.style().xAxisStyle().setLabel("time"); region.style().yAxisStyle().setLabel("Entries"); region.plot(aHisto1D);} {plotter->setCurrentRegionNumber(1); AIDA::IPlotterRegion& region = plotter->currentRegion(); set_region_style(region); region.style().xAxisStyle().setLabel("time"); region.style().yAxisStyle().setLabel("PE"); region.plot(aHisto2D);} plotter->show(); /* FIXME : try to have multiple plotters in the tab stack : {AIDA::IPlotter* plotter = plotterFactory->create("Viewer_2"); if(!plotter) { std::cout << header << "can't get a plotter." << std::endl; return false; } plotter->clearRegions(); plotter->createRegions(1,1,0);} */ plotter->interact(); delete plotter; delete plotterFactory; return true; } ////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////// int main(int aArgc,char** aArgv) ////////////////////////////////////////////////////////////////////////////// //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!// { std::string program = "ELYSE_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 ////////////////////////////////////////// //////////////////////////////////////////////////////// std::string fmt = (aArgc>=2?aArgv[1]:"root"); //file extension AIDA::ITree* tree = treeFactory->create("ELYSE."+fmt,fmt,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; } //FIXME : CINT can't handle dynamic_cast : //AIDA::ITuple* tuple = (AIDA::ITuple*)object->cast("AIDA::ITuple"); AIDA::ITuple* tuple = dynamic_cast(object); if(!tuple) { std::cout << header << "object not an AIDA::ITuple." << std::endl; return 1; } int coln = tuple->columns(); for(int index=0;indexstart(); while(tuple->next() && (irowgetInt(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 tree; delete treeFactory; plot(*aida,*hits_times,*digits_time_pe); delete aida; return 0; }