epic-lfhcal-tbana/OldStructure/CompareDifferentRunsForSinglePhotonSpectra.C

0001 #include <TROOT.h>
0002 #include <TString.h>
0003 #include <TObject.h>
0004 #include <TObjString.h>
0005 #include <TSystem.h>
0006 #include <TChain.h>
0007 #include <TMath.h>
0008 #include <TVector3.h>
0009 #include <iostream>
0010 #include <fstream>
0011 #include <TParticlePDG.h>
0012 #include <TDatabasePDG.h>
0013 #include <TRandom3.h>
0014
0015 #include <TCanvas.h>
0016 #include <TPad.h>
0017 #include <TH1.h>
0018 #include <TH1D.h>
0019 #include <TH1F.h>
0020 #include <TH2.h>
0021 #include <TH3.h>
0022 #include <TFile.h>
0023 #include <TH2D.h>
0024 #include <TH2F.h>
0025 #include <TString.h>
0026 #include <TDatime.h>
0027 #include <TF1.h>
0028 #include <TF2.h>
0029 #include <THStack.h>
0030 #include <TGraph.h>
0031 #include <TStyle.h>
0032 #include <TGraphAsymmErrors.h>
0033 #include <TLine.h>
0034 #include <TLatex.h>
0035 #include <TArrow.h>
0036 #include <TGraphErrors.h>
0037 #include <TGaxis.h>
0038 #include <TLegend.h>
0039 #include <TFrame.h>
0040 #include <TLorentzVector.h>
0041
0042 #include "PlottingHeader.h"
0043 #include "FittingHeader.h"
0044 #include "CommonVariables.h"
0045
0046 const int gMaxLayers    = 16;
0047 const int gMaxSpecies   = 2;
0048 const int gMaxVoltages  = 10;
0049 const int gMaxBEn       = 11;
0050 const int gMaxGainS     = 7;
0051 const int gMaxChannels  = 8;
0052
0053 struct spectraFitDataPoint{
0054     spectraFitDataPoint(): layer(0), channel(0), vop(0), runnr(0), lgSet(0), hgSet(0), assemNr (0) {}
0055     int layer;
0056     int channel;
0057     float vop;
0058     int runnr;
0059     int lgSet;
0060     int hgSet;
0061     int nSPEPeaks;
0062     double avDiffSPEPeaks;
0063     double avDiffSPEPeaksFit;
0064     int assemNr;
0065 } ;
0066 //__________________________________________________________________________________________________________
0067 //_____________________MAIN function !!! ___________________________________________________________________
0068 //__________________________________________________________________________________________________________
0069 void CompareDifferentRunsForSinglePhotonSpectra( TString configFileName   = "configs/BoardA16_RunNumbers",         // list of runs to analyze, filename + runNumber, should be run for one board at the time
0070                                                 TString outputDir         = "Compare/A16/",
0071                                                 Int_t verbosity           = 0,
0072                                                 TString mappingFile       = "configs/mappingA16_SinglePhoton.txt", // board name here
0073                                                 TString runListFileName   = "configs/ORNL_RunNumbers.txt",         // general list of the runs, no board name
0074                                                 Int_t specialData         = 0                                      // specialData: 0- std. TB data, 1 - local ORNL SPE
0075                                   ){
0076     StyleSettingsThesis("pdf");
0077     SetPlotStyle();
0078     Double_t textSizeRel = 0.04;
0079
0080     Color_t colorReadBoard[8] = { kRed+1, kBlue+1, kCyan+1, kMagenta+1, kOrange, kGreen+1, kPink+5, kViolet-9};
0081     Style_t markerReadBoard[8]  = { 20, 21, 33, 34, 29, 39, 40, 46};
0082
0083
0084       // make output directory
0085     TString dateForOutput = ReturnDateStr();
0086     TString outputDirPlots = Form("%s/%s",outputDir.Data(), dateForOutput.Data());
0087     gSystem->Exec("mkdir -p "+outputDir);
0088     gSystem->Exec("mkdir -p "+outputDirPlots);
0089
0090     // ********************************************************************************************************
0091     // read run list and corresponding settings
0092     // ********************************************************************************************************
0093     std::vector<runInfo> allRuns = readRunInfosFromFile(runListFileName, 1, specialData);
0094
0095     // ********************************************************************************************************
0096     // read folder and name from file
0097     // ********************************************************************************************************
0098     std::vector<Int_t> runnumbers;
0099     std::vector<TString> fileNames;
0100     ifstream in;
0101     in.open(configFileName,ios_base::in);
0102     if (!in) {
0103         std::cout << "ERROR: file " << configFileName.Data() << " not found!" << std::endl;
0104         return;
0105     }
0106
0107     for( TString tempLine; tempLine.ReadLine(in, kTRUE); ) {
0108       // check if line should be considered
0109       if (tempLine.BeginsWith("%") || tempLine.BeginsWith("#")){
0110           continue;
0111       }
0112       if (verbosity > 0) std::cout << tempLine.Data() << std::endl;
0113
0114       // Separate the string according to tabulators
0115       TObjArray *tempArr  = tempLine.Tokenize("\t");
0116       if(tempArr->GetEntries()<1){
0117           if (verbosity > 1) std::cout << "nothing to be done" << std::endl;
0118           delete tempArr;
0119           continue;
0120       } else if (tempArr->GetEntries() == 1 ){
0121           // Separate the string according to space
0122           tempArr       = tempLine.Tokenize(" ");
0123           if(tempArr->GetEntries()<1){
0124               if (verbosity > 1) std::cout << "nothing to be done" << std::endl;
0125               delete tempArr;
0126               continue;
0127           } else if (tempArr->GetEntries() == 1  ) {
0128               if (verbosity > 1) std::cout << ((TString)((TObjString*)tempArr->At(0))->GetString()).Data() << " has not been reset, no value given!" << std::endl;
0129               delete tempArr;
0130               continue;
0131           }
0132       }
0133
0134       TString tempname  = ((TString)((TObjString*)tempArr->At(1))->GetString());
0135       Int_t temprun     = ((TString)((TObjString*)tempArr->At(0))->GetString()).Atoi();
0136       runnumbers.push_back(temprun);
0137       fileNames.push_back(tempname);
0138     }
0139
0140     std::cout<<"=============================================================="<<std::endl;
0141     std::vector<runInfo> compRuns;
0142     std::vector<TGraphErrors*>graphsNSPEPeaks;
0143     std::vector<TGraphErrors*>graphsAvDiffSPEPeaks;
0144     std::vector<TGraphErrors*>graphsAvDiffSPEPeaksFit;
0145
0146     Int_t runsToComp = (Int_t)runnumbers.size();
0147
0148     for (Int_t r = 0; r < runsToComp; r++){
0149
0150       Int_t indexCRun = findCurrentRun(allRuns, runnumbers.at(r));
0151       runInfo currentRunInfo;
0152       if (indexCRun < 0){
0153         std::cout << "run not in current list of runs, provided" << std::endl;
0154         return;
0155       } else {
0156         std::cout << "Run "<< runnumbers.at(r) << "\t found at index " << indexCRun << std::endl;
0157         currentRunInfo = GetRunInfoObject(allRuns,indexCRun);
0158       }
0159       compRuns.push_back(currentRunInfo);
0160       if (specialData == 1){
0161         std::cout <<  runnumbers.at(r) << "\t" << fileNames.at(r).Data() << "\t" << compRuns.at(r).vop << "\t"<< compRuns.at(r).assemblyNr << std::endl;
0162       } else {
0163         std::cout <<  runnumbers.at(r) << "\t" << fileNames.at(r).Data() << "\t" << compRuns.at(r).vop  << std::endl;
0164       }
0165
0166       TFile* tempFile = new TFile(fileNames.at(r).Data(), "OPEN");
0167       if (tempFile->IsZombie()){
0168           std::cout << fileNames.at(r).Data() << " is broken, please remove from list or fix!" << std::endl;
0169       }
0170       TGraphErrors* tempGraph = (TGraphErrors*)tempFile->Get("graphnSPEPeaks_HG_channels");
0171       graphsNSPEPeaks.push_back(tempGraph);
0172       tempGraph = (TGraphErrors*)tempFile->Get("graphAvDiffSPEPeaks_HG_channels");
0173       graphsAvDiffSPEPeaks.push_back(tempGraph);
0174       tempGraph = (TGraphErrors*)tempFile->Get("graphAvDiffSPEPeaksFit_HG_channels");
0175       graphsAvDiffSPEPeaksFit.push_back(tempGraph);
0176     }
0177     std::cout<<"=============================================================="<<std::endl;
0178
0179     Double_t rangeDiffPeaks[2] = {0.};
0180     std::vector<spectraFitDataPoint> spectraFitPoints;
0181     for (Int_t r = 0; r < runsToComp; r++){
0182       std::cout <<  runnumbers.at(r) << "\t" << fileNames.at(r).Data() << "\t" << compRuns.at(r).species.Data() << "\t" << compRuns.at(r).vop << std::endl;
0183       // if (graphsNSPEPeaks.at(r)->GetN() !=graphsAvDiffSPEPeaks.at(r)->GetN() || graphsNSPEPeaks.at(r)->GetN() !=graphsAvDiffSPEPeaksFit.at(r)->GetN() ){
0184       //   std::cout << "Graphs don't have same size ... something is wrong!" << std::endl;
0185       //   continue;
0186       // }
0187       for (Int_t e = 0; e < (Int_t)graphsNSPEPeaks.at(r)->GetN(); e++){
0188         spectraFitDataPoint currFit;
0189         currFit.layer   = (int)(graphsNSPEPeaks.at(r)->GetX()[e]/10);
0190         currFit.channel = (int)(graphsNSPEPeaks.at(r)->GetX()[e]-10*currFit.layer);
0191         currFit.vop     = (float)(compRuns.at(r).vop);
0192         currFit.runnr   = (int)(compRuns.at(r).runNr);
0193         currFit.lgSet   = (float)(compRuns.at(r).lgSet);
0194         currFit.hgSet   = (float)(compRuns.at(r).hgSet);
0195         currFit.nSPEPeaks          = (int)(graphsNSPEPeaks.at(r)->GetY()[e]);
0196         currFit.avDiffSPEPeaks     = (double)(graphsAvDiffSPEPeaks.at(r)->GetY()[e]);
0197         currFit.avDiffSPEPeaksFit  = (double)(graphsAvDiffSPEPeaksFit.at(r)->GetY()[e]);
0198         if (specialData == 1) currFit.assemNr = (int)(compRuns.at(r).assemblyNr);
0199 //         std::cout << graphsMPV.at(r)->GetX()[e] << "\t l: " << currMip.layer << "\t c: " << currMip.channel << "\t vop: " << currMip.vop << std::endl;
0200         spectraFitPoints.push_back(currFit);
0201         if (currFit.avDiffSPEPeaks < rangeDiffPeaks[0]) rangeDiffPeaks[0] = currFit.avDiffSPEPeaks;
0202         if (currFit.avDiffSPEPeaks > rangeDiffPeaks[1]) rangeDiffPeaks[1] = currFit.avDiffSPEPeaks;
0203       }
0204     }
0205
0206     Int_t gMaxAssemblies = 100;
0207     TGraphErrors* graphAvSPEDiffVsVop[gMaxAssemblies][gMaxChannels];
0208     TGraphErrors* graphAvSPEDiffFitVsVop[gMaxAssemblies][gMaxChannels];
0209
0210     for (Int_t a = 0; a < gMaxAssemblies; a++){
0211       for (Int_t ch=0; ch<gMaxChannels; ch++){
0212         graphAvSPEDiffVsVop[a][ch]     = new TGraphErrors();
0213         graphAvSPEDiffFitVsVop[a][ch]  = new TGraphErrors();
0214       }
0215     }
0216     for(Int_t f = 0; f<(Int_t)spectraFitPoints.size(); f++){
0217       Int_t ch      = spectraFitPoints.at(f).channel-1;
0218       Int_t nPeaks  = spectraFitPoints.at(f).nSPEPeaks;
0219       Int_t assem   = spectraFitPoints.at(f).layer;
0220       if (specialData == 1){
0221         assem = spectraFitPoints.at(f).assemNr;
0222       }
0223       graphAvSPEDiffVsVop[assem][ch]->SetPoint(graphAvSPEDiffVsVop[assem][ch]->GetN(), spectraFitPoints.at(f).vop, spectraFitPoints.at(f).avDiffSPEPeaks);
0224       graphAvSPEDiffFitVsVop[assem][ch]->SetPoint(graphAvSPEDiffFitVsVop[assem][ch]->GetN(), spectraFitPoints.at(f).vop, spectraFitPoints.at(f).avDiffSPEPeaksFit);
0225     }
0226
0227     for (Int_t a = 0; a < gMaxAssemblies; a++){
0228       for(Int_t ch=0; ch<gMaxChannels; ch++){
0229 //         std::cout << "A" << a <<" channel: " << ch << std::endl;
0230         if( graphAvSPEDiffVsVop[a][ch]->GetN() !=0)      graphAvSPEDiffVsVop[a][ch]->Sort();
0231         if( graphAvSPEDiffFitVsVop[a][ch]->GetN() !=0)   graphAvSPEDiffFitVsVop[a][ch]->Sort();
0232       }
0233     }
0234
0235     // ********************************************************************************************************
0236     // plotting
0237     // ********************************************************************************************************
0238
0239     TCanvas* canvas1DSimple = new TCanvas("canvas1DSimple","",0,0,700,500);  // gives the page size
0240     DefaultCancasSettings( canvas1DSimple, 0.08, 0.02, 0.02, 0.08);
0241
0242
0243     TH1D* histDummyAvDiffSPEvsVop = new TH1D("histDummyAvDiffSPEvsVop","", 55, 1.5, 7 );
0244     SetStyleHistoTH1ForGraphs( histDummyAvDiffSPEvsVop, "#it{V}_{ov} (V)", "#mu(#Delta_{SPE}) (HG ADC)",0.85* textSizeRel, textSizeRel, 0.85* textSizeRel, textSizeRel, 0.9, 0.95);
0245     histDummyAvDiffSPEvsVop->GetYaxis()->SetRangeUser(rangeDiffPeaks[0], rangeDiffPeaks[1]*1.5);
0246     histDummyAvDiffSPEvsVop->Draw();
0247     Bool_t lenteredAll[2]  = {0, 0};
0248     Bool_t plottedA      = 0;
0249     TLegend* legend = GetAndSetLegend2( 0.12, 0.95-3*textSizeRel, 0.7, 0.95,textSizeRel, 4, "RB-Channel:", 42,0.1);
0250
0251     for (Int_t a = 0; a < gMaxAssemblies; a++){
0252       for(Int_t ch=0; ch<gMaxChannels; ch++){
0253         if(graphAvSPEDiffVsVop[a][ch]->GetN() != 0 ) {
0254           SetMarkerDefaultsTGraphErr( graphAvSPEDiffVsVop[a][ch], markerReadBoard[ch], 0.8, colorReadBoard[ch],markerReadBoard[ch]);
0255           graphAvSPEDiffVsVop[a][ch]->Draw("p,e,same");
0256           legend->AddEntry(graphAvSPEDiffVsVop[a][ch], Form("A%02d, %d",a, ch+1), "p");
0257         }
0258       }
0259     }
0260     // DrawLatex(0.95, 0.92, Form("%s", labelTB.Data()), true, textSizeRel, 42);
0261     legend->Draw();
0262     canvas1DSimple->SaveAs(Form("%s/AvDiffSPEPeaks.pdf",outputDirPlots.Data()));
0263
0264     TH1D* histDummyAvDiffSPEvsVopF = new TH1D("histDummyAvDiffSPEvsVopF","", 55, 1.5, 7 );
0265     SetStyleHistoTH1ForGraphs( histDummyAvDiffSPEvsVopF, "#it{V}_{ov} (V)", "#mu(#Delta_{SPE,fit}) (HG ADC)",0.85* textSizeRel, textSizeRel, 0.85* textSizeRel, textSizeRel, 0.9, 0.95);
0266     histDummyAvDiffSPEvsVopF->GetYaxis()->SetRangeUser(rangeDiffPeaks[0], rangeDiffPeaks[1]*1.5);
0267     histDummyAvDiffSPEvsVopF->Draw();
0268
0269     for (Int_t a = 0; a < gMaxAssemblies; a++){
0270       for(Int_t ch=0; ch<gMaxChannels; ch++){
0271         if( graphAvSPEDiffFitVsVop[a][ch]->GetN() != 0 ){
0272           SetMarkerDefaultsTGraphErr( graphAvSPEDiffFitVsVop[a][ch], markerReadBoard[ch], 0.8, colorReadBoard[ch],markerReadBoard[ch]);
0273           graphAvSPEDiffFitVsVop[a][ch]->Draw("p,e,same");
0274         }
0275
0276       }
0277     }
0278     // DrawLatex(0.95, 0.92, Form("%s", labelTB.Data()), true, textSizeRel, 42);
0279     legend->Draw();
0280     canvas1DSimple->SaveAs(Form("%s/AvDiffSPEPeaksFits.pdf",outputDirPlots.Data()));
0281 }