epic-lfhcal-tbana/NewStructure/ComparisonInjection.cc

0001 #include "ComparisonInjection.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 ComparisonInjection::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 ComparisonInjection::ProcessInjectionCompare(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   // ******************************************************************************************
0160   // Iterate over all entries (runs) in the calib tree
0161   // ******************************************************************************************
0162   for(int ientry=0; ientry<entries;ientry++){
0163     TsetupIn->GetEntry(ientry);
0164     TcalibIn->GetEntry(ientry);
0165
0166     // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0167     // set global iterator for runs to first run number in list to obtain beam-line, dates...
0168     if (ientry==0) firstRunNr = calib.GetRunNumber();
0169     itRun = ri.find(calib.GetRunNumber());
0170     double set_rf     = (double)itRun->second.rf;
0171     double set_cf     = (double)itRun->second.cf;
0172     double set_cc     = (double)itRun->second.cc;
0173     double set_cfcomp = (double)itRun->second.cfcomp;
0174     std::cout <<calib.GetRunNumber() << "\t" << set_rf << "\t" << set_cf << "\t" << set_cc << "\t" << set_cfcomp << std::endl;
0175
0176     // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0177     // Set X-values according to option
0178     // Xaxis:
0179     //        0 - Run number dependence
0180     //        1 - Operational Voltage dependence
0181     //        2 - Time dependence
0182     // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0183     Int_t runNumber = calib.GetRunNumber();
0184     Xvalue=calib.GetRunNumber();
0185     if(Xvalue<Xmin) Xmin=Xvalue;
0186     if(Xvalue>Xmax) Xmax=Xvalue;
0187
0188     // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0189     // Initialize calib summary
0190     // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0191     calib.PrintGlobalInfo();
0192     CalibSummary aSum = CalibSummary(nRun, runNumber,calib.GetVop());
0193
0194     // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0195     // Reading additional summary histos from 2nd file
0196     // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0197     TFile* tempFile = nullptr;
0198     if (nRun < (int)RootInputNames.size()){
0199       std::cout << "reading hist file: " <<RootInputNames[nRun].Data() << " expanded list setting: " << expandedList << std::endl;
0200       tempFile      = new TFile(RootInputNames[nRun].Data(),"READ");
0201     }
0202
0203     // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0204     // Loop over all cells in the calib object for trending plots
0205     // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0206
0207     for(itcalib=calib.begin(); itcalib!=calib.end(); ++itcalib){
0208
0209       // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0210       // Reading additional cell histos from 2nd file
0211       // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0212       TProfile* profCellWave= nullptr;
0213       TProfile* profCellTOA= nullptr;
0214       TProfile* profCellTOT= nullptr;
0215       // std::cout << "celld id: " << itcalib->first << "\t "<< calib.GetPedestalMeanH(itcalib->first);
0216       // reading Waveform
0217       profCellWave     = (TProfile*)tempFile->Get(Form("IndividualCells/waveform1DinjCellID%i",itcalib->first));
0218       // reading ToA
0219       profCellTOA      = (TProfile*)tempFile->Get(Form("IndividualCells/TOAinjCellID%i",itcalib->first));
0220       // reading ToT
0221       profCellTOT      = (TProfile*)tempFile->Get(Form("IndividualCells/TOTinjCellID%i",itcalib->first));
0222       // std::cout <<"\t Wave: "<< profCellWave << "\t TOA: " << profCellTOA << "\t TOT " << profCellTOT << std::endl;
0223
0224       // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0225       // fill calib summary object for specific cell
0226       // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0227       aSum.Fill(itcalib->second);
0228
0229       // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0230       // fill trending object for a single cell
0231       // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0232       // check if iterator points to end of map
0233       itrend=trend.find(itcalib->first);
0234       if(itrend!=trend.end()){
0235         // fill injection hists
0236         itrend->second.Fill(Xvalue,itcalib->second, (int)calib.GetRunNumber(), (double)calib.GetVop());
0237         itrend->second.FillInjection(Xvalue,calib.GetPedestalMeanH(itcalib->first),(int)calib.GetRunNumber(),
0238                                      profCellWave, profCellTOA, profCellTOT,
0239                                      set_rf, set_cf, set_cfcomp, set_cc );
0240       // create new TileTrend object if not yet available in map
0241       } else {
0242         TileTrend atrend=TileTrend(itcalib->first,0, 3);
0243         // fill injection hists
0244         atrend.Fill(Xvalue,itcalib->second, (int)calib.GetRunNumber(), (double)calib.GetVop());
0245         atrend.FillInjection(Xvalue,calib.GetPedestalMeanH(itcalib->first),(int)calib.GetRunNumber(),
0246                              profCellWave, profCellTOA, profCellTOT,
0247                              set_rf, set_cf, set_cfcomp, set_cc);
0248         // append TileTrend object to map
0249         trend[itcalib->first]=atrend;
0250       }
0251     } // end loop over cells in the calib object
0252
0253     // append CalibSummary object to map
0254     sumCalibs[nRun]=aSum;
0255     // close additional files opened
0256     if (expandedList){
0257       tempFile->Close();
0258     }
0259     // increase run-counter
0260     nRun++;
0261   } // end loop over entries (runs) in calib tree
0262
0263   // ******************************************************************************************
0264   // Print summary of calib runs
0265   // ******************************************************************************************
0266   std::cout << "Calibs summary: "<< sumCalibs.size() << std::endl;
0267   int globalStatus = 0;
0268   for(isumCalibs=sumCalibs.begin(); isumCalibs!=sumCalibs.end(); ++isumCalibs){
0269     int calibstatus = isumCalibs->second.Analyse(debug);
0270     if (globalStatus < calibstatus) globalStatus = calibstatus;
0271   }
0272   std::cout << "Global calib status: " << globalStatus << std::endl;
0273
0274   if (globalStatus == 0){
0275     std::cout << "!!!!!!!!!!!!!!!!!!!!!  ATTENTION !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl;
0276     std::cout << "Aborting plotting: none of the files has either ped or mip scales filled" << std::endl;
0277   }
0278
0279   // ******************************************************************************************
0280   // Set X axis title and ranges
0281   // ******************************************************************************************
0282   if (Xaxis == 0){
0283     Xmin= Xmin-10;
0284     Xmax= Xmax+10;
0285   }
0286
0287   TString             xaxisTitle = "";
0288   if (Xaxis==0)       xaxisTitle = "Run Nr. ";
0289
0290   for(itrend=trend.begin(); itrend!=trend.end(); ++itrend){
0291     // sort graphs
0292     itrend->second.Sort();
0293     // set x axis title for trending graphs
0294     itrend->second.SetXAxisTitle(xaxisTitle);
0295     // write graphs for each cell to output
0296     itrend->second.Write(RootOutput);
0297   }
0298
0299   itRun = ri.find(firstRunNr);
0300
0301   bool isSameVoltage    = true;
0302   double commonVoltage  = 0;
0303   bool isSameRF         = true;
0304   bool isSameCF         = true;
0305   bool isSameCFcomp     = true;
0306   bool isSameCC         = true;
0307   double commonRF       = 0;
0308   double commonCF       = 0;
0309   double commonCFcomp   = 0;
0310   double commonCC       = 0;
0311   itrend=trend.begin();
0312   for (int rc = 0; rc < itrend->second.GetNRuns() && rc < 30; rc++ ){
0313     if (rc == 0){
0314       commonVoltage = itrend->second.GetVoltage(rc);
0315       commonRF      = itrend->second.GetRF(rc);
0316       commonCF      = itrend->second.GetCF(rc);
0317       commonCC      = itrend->second.GetCC(rc);
0318       commonCFcomp  = itrend->second.GetCFComp(rc);
0319     } else {
0320       if (commonVoltage != itrend->second.GetVoltage(rc))   isSameVoltage = false;
0321       if (commonRF != itrend->second.GetRF(rc))        isSameRF = false;
0322       if (commonCF != itrend->second.GetCF(rc))        isSameCF = false;
0323       if (commonCFcomp != itrend->second.GetCFComp(rc))    isSameCFcomp = false;
0324       if (commonCC != itrend->second.GetCC(rc))        isSameCC = false;
0325     }
0326     std::cout << itrend->second.GetRunNr(rc) << "\t" <<  itrend->second.GetRF(rc) << "\t" << itrend->second.GetCF(rc) << "\t" << itrend->second.GetCC(rc) << "\t" << itrend->second.GetCFComp(rc) << "\tRF calc:"<< ReturnRFValue(itrend->second.GetRF(rc))<< std::endl;
0327
0328   }
0329   std::cout << "Show all common settings" << std::endl;
0330   std::cout << "Common Voltage: "<< isSameVoltage << "\t"  << commonVoltage
0331             << "\t RF: " << isSameRF << "\t"  << commonRF
0332             << "\t CF: " << isSameCF << "\t"  << commonCF
0333             << "\t CFcomp: " << isSameCFcomp << "\t"  << commonCFcomp
0334             << "\t CC: " << isSameCC << "\t"  << commonCC
0335             << std::endl;
0336
0337   RunInfo commonRunInfo;
0338   commonRunInfo.nFPGA       = itRun->second.nFPGA;
0339   commonRunInfo.nASIC       = itRun->second.nASIC;
0340   commonRunInfo.energy      = itRun->second.energy;
0341   commonRunInfo.pdg         = itRun->second.pdg;
0342   commonRunInfo.species     = itRun->second.species;
0343   commonRunInfo.beamline    = itRun->second.beamline;
0344   commonRunInfo.facility    = itRun->second.facility;
0345   commonRunInfo.readout     = itRun->second.readout;
0346   commonRunInfo.month       = itRun->second.month;
0347   commonRunInfo.year        = itRun->second.year;
0348   commonRunInfo.samples     = itRun->second.samples;
0349   commonRunInfo.samples     = itRun->second.vbr;
0350   if (isSameRF)
0351     commonRunInfo.rf        = commonRF;
0352   else
0353     commonRunInfo.rf        = -10.;
0354   if (isSameCF)
0355     commonRunInfo.cf        = commonCF;
0356   else
0357     commonRunInfo.cf        = -10.;
0358   if (isSameCFcomp)
0359     commonRunInfo.cfcomp    = commonCFcomp;
0360   else
0361     commonRunInfo.cfcomp    = -10.;
0362   if (isSameCC)
0363     commonRunInfo.cc        = commonCC;
0364   else
0365     commonRunInfo.cc        = -10.;
0366   if (isSameVoltage)
0367     commonRunInfo.vop       = commonVoltage;
0368   else
0369     commonRunInfo.vop       = -10.;
0370   //******************************************************************************
0371   // plotting overview for each run overlayed
0372   //******************************************************************************
0373   Int_t textSizePixel   = 30;
0374   Float_t textSizeRel   = 0.04;
0375   TCanvas* canvas1DRunsOverlay = new TCanvas("canvas1DRunsOverlay","",0,0,1450,1300);  // gives the page size
0376   DefaultCanvasSettings( canvas1DRunsOverlay, 0.075, 0.015, 0.025, 0.09);
0377
0378   PlotCalibRunOverlay( canvas1DRunsOverlay, 0, sumCalibs, textSizeRel,
0379                       Form("%s/HGPedSummary_RunOverlay.%s",OutputNameDirPlots.Data(),plotSuffix.Data()), commonRunInfo,"", debug);
0380   PlotCalibRunOverlay( canvas1DRunsOverlay, 1, sumCalibs, textSizeRel,
0381                       Form("%s/HGPedWidthSummary_RunOverlay.%s",OutputNameDirPlots.Data(),plotSuffix.Data()), commonRunInfo,"", debug);
0382
0383
0384   // plotting individual layers/asics
0385   DetConf::Type detConf = DetConf::Type::Asic;
0386   // DetConf::Type detConf = DetConf::Type::Single8M;
0387   // DetConf::Type detConf = DetConf::Type::Dual8M;
0388   MultiCanvas panelPlot2D(detConf, "Injection");
0389   bool init2D = panelPlot2D.Initialize(2);
0390
0391   // panelPlot2D.PlotTrending(trend, 0, Xmin,Xmax, OutputNameDirPlots, "PedADC", plotSuffix, commonRunInfo, ExtPlot );
0392   panelPlot2D.PlotRunOverlayProfile(trend, nRun, 1, 0, 18*7, -10, -10000, OutputNameDirPlots, "WaveOverlay", plotSuffix, commonRunInfo, ExtPlot );
0393   // panelPlot2D.PlotRunOverlayProfile(trend, nRun, 3, 0, 18*7, 0,-10000, OutputNameDirPlots, "TOTOverlay", plotSuffix, commonRunInfo, ExtPlot );
0394   // panelPlot2D.PlotRunOverlayProfile(trend, nRun, 2, 0, 18*7, 0,1024, OutputNameDirPlots, "TOAOverlay", plotSuffix, commonRunInfo, ExtPlot );
0395
0396   detConf = DetConf::Type::SingleTile;
0397   MultiCanvas panelSingleTile(detConf, "InjectionTile");
0398   bool initSngle = panelSingleTile.Initialize(1);
0399
0400   panelSingleTile.PlotRunOverlayProfile(trend, nRun, 1, 0, 18*7, -10, -10000, OutputNameDirPlots, "TileWaveOverlay", plotSuffix, commonRunInfo, ExtPlot );
0401   return status;
0402 }
0403
0404
0405 // ===========================================================================================
0406 // Create the output file
0407 // ===========================================================================================
0408 bool ComparisonInjection::CreateOutputRootFile(void){
0409   if(Overwrite){
0410     RootOutput=new TFile(RootOutputName.Data(),"RECREATE");
0411   } else{
0412     RootOutput = new TFile(RootOutputName.Data(),"CREATE");
0413   }
0414   if(RootOutput->IsZombie()){
0415     std::cout<<"Error opening '"<<RootOutput<<"'no reachable path? Exist without force mode to overwrite?..."<<std::endl;
0416     return false;
0417   }
0418   return true;
0419 }