epic-lfhcal-tbana/NewStructure/ComparisonWaveform.cc

0001 #include "ComparisonWaveform.h"
0002 #include <vector>
0003 #include "TROOT.h"
0004 #ifdef __APPLE__
0005 #include <unistd.h>
0006 #endif
0007 #include "TF1.h"
0008 #include "TFitResult.h"
0009 #include "TFitResultPtr.h"
0010 #include "TH1D.h"
0011 #include "TH2D.h"
0012 #include "TProfile.h"
0013 #include "TChain.h"
0014 #include "TileSpectra.h"
0015 #include "TileTrend.h"
0016 #include "CalibSummary.h"
0017 #include "MultiCanvas.h"
0018 #include "CommonHelperFunctions.h"
0019 #include "PlotHelper.h"
0020
0021
0022 // ===========================================================================================
0023 // Check input files and global settings
0024 // ===========================================================================================
0025 bool ComparisonWaveform::CheckAndOpenIO(void){
0026
0027   int matchingbranch;
0028
0029   // *****************************************************************************************
0030   // Reading files from a text file
0031   // *****************************************************************************************
0032   if(!InputListName.IsNull()){
0033     // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0034     // text file with 2 files per line 1 full file & 1 histo file
0035     // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0036     std::cout << "You need to provide data tree file with the setup & calib from the injection runs & a histo file" << std::endl;
0037     std::fstream dummyTxt;
0038     dummyTxt.open(InputListName.Data(),std::ios::in);
0039     if(!dummyTxt.is_open()){
0040       std::cout<<"Error opening "<<InputListName.Data()<<", does the file exist?"<<std::endl;
0041     }
0042     std::string dummyRootCalibName;
0043     std::string dummyRootHistName;
0044     // set first root file names
0045     dummyTxt>>dummyRootCalibName >> dummyRootHistName;
0046
0047     int goodcalib;
0048     int goodsetup;
0049     while(dummyTxt.good()){
0050       std::cout << dummyRootCalibName.data() << "\t" << dummyRootHistName.data() << std::endl;
0051
0052       // check that files exist and can be opened
0053       TFile dummyRootCalib=TFile(dummyRootCalibName.c_str(),"READ");
0054       if(dummyRootCalib.IsZombie()){
0055         std::cout<<"Error opening '"<<dummyRootCalibName<<", does the file exist?"<<std::endl;
0056         return false;
0057       }
0058       dummyRootCalib.Close();
0059       TFile dummyRootHist=TFile(dummyRootHistName.c_str(),"READ");
0060       if(dummyRootHist.IsZombie()){
0061         std::cout<<"Error opening '"<<dummyRootHistName<<", does the file exist?"<<std::endl;
0062         return false;
0063       }
0064       dummyRootHist.Close();
0065
0066       // Add file-name to setup and calib chain as well string-vector
0067       AddInputFile(dummyRootHistName);
0068       goodsetup=TsetupIn->AddFile(dummyRootCalibName.c_str());
0069       goodcalib=TcalibIn->AddFile(dummyRootCalibName.c_str());
0070       if(goodcalib==0){
0071           std::cout<<"Issues retrieving Calib tree from "<<dummyRootCalibName<<", file is ignored"<<std::endl;
0072         }
0073       if(goodsetup==0){
0074         std::cout<<"Issues retrieving Setup tree from "<<dummyRootCalibName<<", file is ignored"<<std::endl;
0075       }
0076       // set next root file names
0077       dummyTxt>>dummyRootCalibName >> dummyRootHistName;
0078     }
0079   }
0080   // *****************************************************************************************
0081   // Setup Output files
0082   // *****************************************************************************************
0083   if(RootOutputName.IsNull()){
0084     return false;
0085   } else {
0086     if(!CreateOutputRootFile()){
0087       return false;
0088     }
0089   }
0090
0091   // *****************************************************************************************
0092   // Setup TChain of setups and calibrations
0093   // *****************************************************************************************
0094   // intialize global variable setup
0095   setup=Setup::GetInstance();
0096   std::cout<<"Setup add "<<setup<<std::endl;
0097   matchingbranch=TsetupIn->SetBranchAddress("setup",&rswptr);
0098   if(matchingbranch<0){
0099     std::cout<<"Error retrieving Setup info from the tree"<<std::endl;
0100     return false;
0101   }
0102   std::cout<<"Entries "<<TsetupIn->GetEntries()<<std::endl;
0103   TsetupIn->GetEntry(0);
0104   setup->Initialize(*rswptr);
0105   // initialize calib with the correct branch
0106   matchingbranch=TcalibIn->SetBranchAddress("calib",&calibptr);
0107   if(matchingbranch<0){
0108     std::cout<<"Error retrieving calibration info from the tree"<<std::endl;
0109     return false;
0110   }
0111
0112   return true;
0113 }
0114
0115 // ===========================================================================================
0116 // Main function of this calibration comparison
0117 // ===========================================================================================
0118 bool ComparisonWaveform::ProcessWaveformCompare(void){
0119   // *****************************************************************************************
0120   // plotting settings
0121   // *****************************************************************************************
0122   gSystem->Exec("mkdir -p "+OutputNameDirPlots);
0123   if (ExtPlot > 0) gSystem->Exec("mkdir -p "+OutputNameDirPlots+"/SingleLayer");
0124   StyleSettingsBasics("pdf");
0125   SetPlotStyle();
0126
0127   // *****************************************************************************************
0128   // Some general setup
0129   // *****************************************************************************************
0130   bool status=true;
0131   // enbale implitcit root multithreading
0132   ROOT::EnableImplicitMT();
0133   // get nuber of entires from Calib tree (how many runs do we have)
0134   int entries=TcalibIn->GetEntries();
0135   std::cout << "Entries in calib tree: " << entries << std::endl;
0136
0137   // *****************************************************************************************
0138   // global variable setup, common iterators and ranges
0139   // ******************************************************************************************
0140   std::map<int, TileTrend> trend;
0141   std::map<int, TileTrend>::iterator itrend;
0142   std::map<int, TileCalib>::const_iterator itcalib;
0143
0144   std::map<int, CalibSummary> sumCalibs;
0145   std::map<int, CalibSummary>::iterator isumCalibs;
0146
0147   double Xvalue;
0148   double Xmin= 9999.;
0149   double Xmax=-9999.;
0150   int nRun = 0;
0151
0152   // ******************************************************************************************
0153   // ************* Get run data base to potentially obtain more information from file *********
0154   // ******************************************************************************************
0155   std::map<int,RunInfo> ri=readRunInfosFromFile(RunListInputName.Data(),debug,0);
0156   std::map<int,RunInfo>::iterator itRun; // basic infos
0157   int firstRunNr    = -1;
0158
0159   std::vector<RunInfo> runList;
0160
0161   // ******************************************************************************************
0162   // Iterate over all entries (runs) in the calib tree
0163   // ******************************************************************************************
0164   for(int ientry=0; ientry<entries;ientry++){
0165     TsetupIn->GetEntry(ientry);
0166     TcalibIn->GetEntry(ientry);
0167
0168     // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0169     // set global iterator for runs to first run number in list to obtain beam-line, dates...
0170     if (ientry==0) firstRunNr = calib.GetRunNumber();
0171     itRun = ri.find(calib.GetRunNumber());
0172     runList.push_back(itRun->second);
0173     double set_rf     = (double)itRun->second.rf;
0174     double set_cf     = (double)itRun->second.cf;
0175     double set_cc     = (double)itRun->second.cc;
0176     double set_cfcomp = (double)itRun->second.cfcomp;
0177     double set_injec  = -10000.;
0178     std::cout <<calib.GetRunNumber() << "\t" << set_rf << "\t" << set_cf << "\t" << set_cc << "\t" << set_cfcomp << "\t" << set_injec<< std::endl;
0179
0180     // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0181     // Set X-values according to option
0182     // Xaxis:
0183     //        0 - Run number dependence
0184     //        1 - Operational Voltage dependence
0185     //        2 - Time dependence
0186     // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0187     Int_t runNumber = calib.GetRunNumber();
0188     Xvalue=calib.GetRunNumber();
0189     if(Xvalue<Xmin) Xmin=Xvalue;
0190     if(Xvalue>Xmax) Xmax=Xvalue;
0191
0192     // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0193     // Initialize calib summary
0194     // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0195     calib.PrintGlobalInfo();
0196     CalibSummary aSum = CalibSummary(nRun, runNumber,calib.GetVop(), 0);
0197
0198     // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0199     // Reading additional summary histos from 2nd file
0200     // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0201     TFile* tempFile = nullptr;
0202     if (nRun < (int)RootInputNames.size()){
0203       std::cout << "reading hist file: " <<RootInputNames[nRun].Data() << " expanded list setting: " << expandedList << std::endl;
0204       tempFile      = new TFile(RootInputNames[nRun].Data(),"READ");
0205     }
0206
0207     // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0208     // Loop over all cells in the calib object for trending plots
0209     // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0210
0211     for(itcalib=calib.begin(); itcalib!=calib.end(); ++itcalib){
0212
0213       // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0214       // Reading additional cell histos from 2nd file
0215       // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0216       TProfile* profCellWave= nullptr;
0217       TH1D* histCellTOA= nullptr;
0218       TH1D* histCellTOT= nullptr;
0219
0220       // std::cout << "celld id: " << itcalib->first << "\t "<< calib.GetPedestalMeanH(itcalib->first);
0221       // reading Waveform
0222       profCellWave     = (TProfile*)tempFile->Get(Form("IndividualCells/waveform1DfullCellID%i",itcalib->first));
0223       // reading ToA
0224       histCellTOA      = (TH1D*)tempFile->Get(Form("IndividualCells/hspectrafullTOACellID%i",itcalib->first));
0225       // reading ToT
0226       histCellTOT      = (TH1D*)tempFile->Get(Form("IndividualCells/hspectrafullTOTCellID%i",itcalib->first));
0227       // std::cout <<"\t Wave: "<< profCellWave << "\t TOA: " << histCellTOA << "\t TOT " << histCellTOT << std::endl;
0228       int entries = 0;
0229       if (profCellWave)
0230         entries = profCellWave->GetEntries();
0231
0232       if (entries == 0) continue;
0233       // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0234       // fill calib summary object for specific cell
0235       // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0236       aSum.Fill(itcalib->second);
0237
0238       // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0239       // fill trending object for a single cell
0240       // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0241       // check if iterator points to end of map
0242       itrend=trend.find(itcalib->first);
0243       if(itrend!=trend.end()){
0244         // fill injection hists
0245         itrend->second.Fill(Xvalue,itcalib->second, (int)calib.GetRunNumber(), (double)calib.GetVop(),itRun->second.pdg, -10000., -10000.,itRun->second.energy,itRun->second.temp);
0246         itrend->second.FillExtended(Xvalue,profCellWave->GetEntries(), (int)calib.GetRunNumber(), histCellTOT, histCellTOA, nullptr, profCellWave);
0247         itrend->second.FillHGCROCSetting (set_rf, set_cf, set_cfcomp, set_cc, set_injec);
0248       // create new TileTrend object if not yet available in map
0249       } else {
0250         TileTrend atrend=TileTrend(itcalib->first,0, 2);
0251         // fill injection hists
0252         atrend.Fill(Xvalue,itcalib->second, (int)calib.GetRunNumber(), (double)calib.GetVop(), itRun->second.pdg, -10000., -10000.,itRun->second.energy,itRun->second.temp );
0253         atrend.FillExtended(Xvalue,profCellWave->GetEntries(), (int)calib.GetRunNumber(),histCellTOT,  histCellTOA, nullptr, profCellWave);
0254         atrend.FillHGCROCSetting (set_rf, set_cf, set_cfcomp, set_cc, set_injec);
0255         // append TileTrend object to map
0256         trend[itcalib->first]=atrend;
0257       }
0258     } // end loop over cells in the calib object
0259
0260     // append CalibSummary object to map
0261     sumCalibs[nRun]=aSum;
0262     // close additional files opened
0263     if (expandedList){
0264       tempFile->Close();
0265     }
0266     // increase run-counter
0267     nRun++;
0268   } // end loop over entries (runs) in calib tree
0269
0270   // ******************************************************************************************
0271   // Print summary of calib runs
0272   // ******************************************************************************************
0273   std::cout << "Calibs summary: "<< sumCalibs.size() << std::endl;
0274   int globalStatus = 0;
0275   for(isumCalibs=sumCalibs.begin(); isumCalibs!=sumCalibs.end(); ++isumCalibs){
0276     int calibstatus = isumCalibs->second.Analyse(debug);
0277     if (globalStatus < calibstatus) globalStatus = calibstatus;
0278   }
0279   std::cout << "Global calib status: " << globalStatus << std::endl;
0280
0281   if (globalStatus == 0){
0282     std::cout << "!!!!!!!!!!!!!!!!!!!!!  ATTENTION !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl;
0283     std::cout << "Aborting plotting: none of the files has either ped or mip scales filled" << std::endl;
0284   }
0285
0286   // ******************************************************************************************
0287   // Set X axis title and ranges
0288   // ******************************************************************************************
0289   if (Xaxis == 0){
0290     Xmin= Xmin-10;
0291     Xmax= Xmax+10;
0292   }
0293
0294   TString             xaxisTitle = "";
0295   if (Xaxis==0)       xaxisTitle = "Run Nr. ";
0296
0297   for(itrend=trend.begin(); itrend!=trend.end(); ++itrend){
0298     // sort graphs
0299     itrend->second.Sort();
0300     // set x axis title for trending graphs
0301     itrend->second.SetXAxisTitle(xaxisTitle);
0302     // write graphs for each cell to output
0303     itrend->second.Write(RootOutput);
0304   }
0305
0306   RunInfo commonRunInfo = GetCommonRunInfoFromList(runList);
0307   PrintSettingsRunInfo(commonRunInfo);
0308
0309   //******************************************************************************
0310   // plotting overview for each run overlayed
0311   //******************************************************************************
0312   Int_t textSizePixel   = 30;
0313   Float_t textSizeRel   = 0.04;
0314   TCanvas* canvas1DRunsOverlay = new TCanvas("canvas1DRunsOverlay","",0,0,1450,1300);  // gives the page size
0315   DefaultCanvasSettings( canvas1DRunsOverlay, 0.075, 0.015, 0.025, 0.09);
0316
0317   PlotCalibRunOverlay( canvas1DRunsOverlay, 0, sumCalibs, textSizeRel,
0318                       Form("%s/HGPedSummary_RunOverlay.%s",OutputNameDirPlots.Data(),plotSuffix.Data()), commonRunInfo,"", debug);
0319   PlotCalibRunOverlay( canvas1DRunsOverlay, 1, sumCalibs, textSizeRel,
0320                       Form("%s/HGPedWidthSummary_RunOverlay.%s",OutputNameDirPlots.Data(),plotSuffix.Data()), commonRunInfo,"", debug);
0321
0322   // plotting individual layers/asics
0323   DetConf::Type detConf = DetConf::Type::Asic;
0324   // DetConf::Type detConf = DetConf::Type::Single8M;
0325   // DetConf::Type detConf = DetConf::Type::Dual8M;
0326   MultiCanvas panelPlot2D(detConf, "Waveform");
0327   bool init2D = panelPlot2D.Initialize(2);
0328
0329   // panelPlot2D.PlotTrending(trend, 0, Xmin,Xmax, OutputNameDirPlots, "PedADC", plotSuffix, commonRunInfo, ExtPlot );
0330   panelPlot2D.PlotRunOverlayProfile(trend, nRun, 1, -25, 25*commonRunInfo.samples, -10, -10000, OutputNameDirPlots, "WaveOverlay", plotSuffix, commonRunInfo, ExtPlot );
0331   // panelPlot2D.PlotRunOverlayProfile(trend, nRun, 3, 0, 18*7, 0,-10000, OutputNameDirPlots, "TOTOverlay", plotSuffix, commonRunInfo, ExtPlot );
0332   // panelPlot2D.PlotRunOverlayProfile(trend, nRun, 2, 0, 18*7, 0,1024, OutputNameDirPlots, "TOAOverlay", plotSuffix, commonRunInfo, ExtPlot );
0333
0334   std::vector <int> cellVec;
0335   if (cellList.CompareTo("")!= 0){
0336     std::cout << "cell List set: "  << cellList.Data() << std::endl;
0337     std::fstream cellTxt;
0338     cellTxt.open(cellList.Data(),std::ios::in);
0339     if(!cellTxt.is_open()){
0340       std::cout<<"Error opening "<<cellList.Data()<<", does the file exist?"<<std::endl;
0341     }
0342     while(cellTxt.good()){
0343       TString dummyCell;
0344       // set first root file names
0345       cellTxt>>dummyCell;
0346       std::cout << "\t" << dummyCell.Data() << std::endl;
0347       if (dummyCell.CompareTo("") != 0)
0348         cellVec.push_back(dummyCell.Atoi());
0349     }
0350     std::cout << "registered: " << cellVec.size() << " single cells to be plotted" << std::endl;
0351     for (int i = 0; i < cellVec.size(); i++){
0352         std::cout << cellVec.at(i) << "," ;
0353     }
0354     std::cout << std::endl;
0355   }
0356
0357   if (cellVec.size() > 0){
0358
0359     detConf = DetConf::Type::SingleTile;
0360     MultiCanvas panelSingleTile(detConf, "WaveformTile");
0361     panelSingleTile.SetCellVector(cellVec);
0362     bool initSngle = panelSingleTile.Initialize(1);
0363     panelSingleTile.SetLabelOpt(1);
0364
0365     panelSingleTile.PlotRunOverlayProfile(trend, nRun, 1, -25, 25*commonRunInfo.samples, -10, 1324, OutputNameDirPlots, "TileWaveOverlay", plotSuffix, commonRunInfo, ExtPlot);
0366     panelSingleTile.PlotRunOverlaySpectra(trend, nRun, 0, 0, 4148, OutputNameDirPlots, "TileTOTSpectra", plotSuffix, commonRunInfo, ExtPlot, debug, 0);
0367     panelSingleTile.PlotRunOverlaySpectra(trend, nRun, 1, 0, 1024, OutputNameDirPlots, "TileTOASpectra", plotSuffix, commonRunInfo, ExtPlot, debug, 0);
0368
0369   }
0370   return status;
0371 }
0372
0373
0374 // ===========================================================================================
0375 // Create the output file
0376 // ===========================================================================================
0377 bool ComparisonWaveform::CreateOutputRootFile(void){
0378   if(Overwrite){
0379     RootOutput=new TFile(RootOutputName.Data(),"RECREATE");
0380   } else{
0381     RootOutput = new TFile(RootOutputName.Data(),"CREATE");
0382   }
0383   if(RootOutput->IsZombie()){
0384     std::cout<<"Error opening '"<<RootOutput<<"'no reachable path? Exist without force mode to overwrite?..."<<std::endl;
0385     return false;
0386   }
0387   return true;
0388 }