Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-01-18 09:15:50 

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 #include <TSpectrum.h>
0042 #include <TVirtualFitter.h>
0043 
0044 #include "PlottingHeader.h"
0045 #include "FittingHeader.h"
0046 
0047 const int gMaxChannels = 64;
0048 const int gMaxBoard    = 1;
0049 const int gMaxLayers   = 16;
0050 // Define tree variables
0051 Long64_t gTrID;
0052 Double_t gTRtimeStamp;
0053 Long64_t* gBoard          = new Long64_t[gMaxChannels];
0054 Long64_t* gChannel        = new Long64_t[gMaxChannels];
0055 Long64_t* gLG             = new Long64_t[gMaxChannels];
0056 Long64_t* gHG             = new Long64_t[gMaxChannels];
0057 
0058 // Definition of color/marker arrays for plotting
0059 Color_t colorLayer[16]    = { kRed-7, kBlue+1, kCyan-3, kMagenta+1, kPink-1, kViolet-9, kOrange, kGreen+1,  
0060                               kRed+1, kBlue-9, kCyan+1, kMagenta-7, kOrange+7, kSpring+5, kPink+5, kViolet+8 };
0061 Color_t colorReadBoard[8] = { kRed+1, kBlue+1, kCyan+1, kMagenta+1, kOrange, kGreen+1, kPink+5, kViolet-9};
0062 Style_t markerLayer[16]     = { 20, 25, 33, 28, 29, 37, 41, 46,
0063                                 24, 21, 27, 34, 30, 39, 40, 47};
0064 Style_t markerReadBoard[8]  = { 20, 21, 33, 34, 29, 39, 40, 46};                             
0065 
0066 #include "makeSimplePlotsFromJanusTree.h"
0067 
0068 
0069 
0070 //__________________________________________________________________________________________________________
0071 //_____________________MAIN function !!! ___________________________________________________________________
0072 //__________________________________________________________________________________________________________
0073 void makeSinglePhotonSpectraFitsFromJanusTree(  TString fileName     = "", 
0074                                                 TString outputDir    = "ProcessedData/",
0075                                                 Int_t runnumber      = -1,
0076                                                 ULong_t minNEvent    = 0,
0077                                                 ULong_t maxNEvent    = 0,
0078                                                 Int_t verbosity      = 0,
0079                                                 Int_t dataType       = 0,
0080                                                 TString mappingFile  = "",
0081                                                 Bool_t bDetPlot      = kTRUE,
0082                                                 TString runListFileName = "configs/SPS_RunNumbers.txt"
0083                                               ){
0084                                
0085     //*******************************************************
0086     //********************Running modes:*********************
0087     // dataType = 0       // local testing and parsing
0088     // dataType = 1       // PS test beam 2023
0089     // dataType = 2       // local test setups in small stack
0090   
0091     StyleSettingsThesis("pdf");
0092     SetPlotStyle();
0093     Double_t textSizeRel = 0.04;
0094     
0095       // make output directory
0096     TString dateForOutput = ReturnDateStr();
0097     if (runnumber > -1)
0098       outputDir = Form("%s/Run%05d",outputDir.Data(), runnumber );
0099     TString outputDirPlots = Form("%s/Plots",outputDir.Data());
0100     TString outputDirPlotsDet = Form("%s/Plots/Detailed",outputDir.Data());
0101     gSystem->Exec("mkdir -p "+outputDir);
0102     gSystem->Exec("mkdir -p "+outputDirPlots);
0103     gSystem->Exec("mkdir -p "+outputDirPlotsDet);
0104 
0105     // load tree
0106     TChain *const tt_event = new TChain("tree");
0107     if (fileName.EndsWith(".root")) {                     // are we loading a single root tree?
0108         std::cout << "loading a single root file" << std::endl;
0109         tt_event->AddFile(fileName);
0110         TFile testFile(fileName.Data());
0111         if (testFile.IsZombie()){
0112           std::cout << Form("The file %s is not a root file or doesn't exit, please fix the file path", fileName.Data()) << std::endl;
0113           return;
0114         }
0115         
0116     } else {
0117         std::cout << "please try again this isn't a root file" << std::endl;
0118     } 
0119     if(!tt_event){ std::cout << "tree not found... returning!"<< std::endl; return;}
0120 
0121     // load all branches (see header)
0122     SetBranchAddressesTree(tt_event);
0123 
0124     
0125     Long64_t nEntriesTree                 = tt_event->GetEntries();
0126     std::cout << "Number of events in tree: " << nEntriesTree << std::endl;
0127     if(maxNEvent>0 && maxNEvent<(ULong_t)nEntriesTree){
0128       nEntriesTree = maxNEvent;
0129       std::cout << "Will only analyze first " << maxNEvent << " events in the tree..." << std::endl;
0130     }
0131 
0132     // ********************************************************************************************************
0133     // read run list and corresponding settings
0134     // ********************************************************************************************************
0135     std::vector<runInfo> runs = readRunInfosFromFile(runListFileName, 0 );
0136     Int_t indexCRun = findCurrentRun(runs, runnumber);
0137     runInfo currentRunInfo;
0138     if (indexCRun < 0){
0139       std::cout << "run not in current list of runs, provided" << std::endl;
0140       return;
0141     } else {
0142       std::cout << "Run "<< runnumber << "\t found at index " << indexCRun << std::endl;
0143       currentRunInfo = GetRunInfoObject(runs,indexCRun);
0144       std::cout << "Run " << currentRunInfo.runNr << "\t species: " << currentRunInfo.species << "\t energy: "  << currentRunInfo.energy << "\t Vop: " << currentRunInfo.vop << std::endl;
0145     }
0146     
0147     // ********************************************************************************************************
0148     // Read mapping
0149     //  * association of CAEN board channels with actuall layers in the detectors
0150     //  * CAEN board Nr. 0-2
0151     //  * CAEN board channels 0-64
0152     //  * Layers during TB (0-14)
0153     //  * Readout board channel (1-8)
0154     //  * Readout board column (0-3)
0155     //  * Readout board row (0-1)
0156     //  * Assembly Nr TB (0-19)
0157     // ********************************************************************************************************
0158     Int_t mapping[128][4]                     = {{-1}};      // linear mapping 
0159     Int_t backwardMapping[9][gMaxLayers]      = {{-1}};     // backward associations channels
0160     Int_t backwardMappingBoard[9][gMaxLayers] = {{-1}};     // backward associations board
0161     Bool_t lActive[gMaxLayers]                = {0};        // boolean to check whether layer is active
0162     Bool_t lActiveCh[9]                       = {0};        // boolean to check whether channel on board is active
0163     for (Int_t l = 0; l < gMaxLayers; l++){
0164       for (Int_t c = 0; c < 9 ; c++){
0165         backwardMapping[c][l]       = -1;
0166         backwardMappingBoard[c][l]  = -1;
0167         lActiveCh[c]                = kFALSE;
0168       }
0169       lActive[l] = kFALSE;
0170     }
0171     Int_t nChmapped = 0;
0172     cout << "INFO: You have chosen the given the following config file: " << mappingFile.Data() << endl;
0173     ifstream fileMapping;
0174     fileMapping.open(mappingFile,ios_base::in);
0175     if (!fileMapping) {
0176         cout << "ERROR: settings " << mappingFile.Data() << " not found!" << endl;
0177         return;
0178     }
0179 
0180     // read settings from file
0181     for( TString tempLine; tempLine.ReadLine(fileMapping, kTRUE); ) {
0182         // check if line should be considered
0183         if (tempLine.BeginsWith("%") || tempLine.BeginsWith("#")){
0184             continue;
0185         }
0186         if (verbosity > 0) cout << tempLine.Data() << endl;
0187 
0188         // Separate the string according to tabulators
0189         TObjArray *tempArr  = tempLine.Tokenize("\t");
0190         if(tempArr->GetEntries()<1){
0191             cout << "nothing to be done" << endl;
0192             delete tempArr;
0193             continue;
0194         } else if (tempArr->GetEntries() == 1 ){
0195             // Separate the string according to space
0196             tempArr       = tempLine.Tokenize(" ");
0197             if(tempArr->GetEntries()<1){
0198                 cout << "nothing to be done" << endl;
0199                 delete tempArr;
0200                 continue;
0201             } else if (tempArr->GetEntries() == 1 ) {
0202                 cout << ((TString)((TObjString*)tempArr->At(0))->GetString()).Data() << " has not be reset, no value given!" << endl;
0203                 delete tempArr;
0204                 continue;
0205             }
0206         }
0207         Int_t chCAEN    = ((TString)((TObjString*)tempArr->At(0))->GetString()).Atoi()*64 + ((TString)((TObjString*)tempArr->At(1))->GetString()).Atoi();
0208         Int_t layerMod  = ((TString)((TObjString*)tempArr->At(2))->GetString()).Atoi();
0209         Int_t assemblyMod  = ((TString)((TObjString*)tempArr->At(3))->GetString()).Atoi();
0210         Int_t chBoard   = ((TString)((TObjString*)tempArr->At(4))->GetString()).Atoi();
0211         Int_t rowBoard  = ((TString)((TObjString*)tempArr->At(5))->GetString()).Atoi();
0212         Int_t colBoard  = ((TString)((TObjString*)tempArr->At(6))->GetString()).Atoi();
0213         
0214         if (verbosity > 0) std::cout << "-->" << chCAEN << "\t" << layerMod << "\t"<< chBoard << "\t" << rowBoard << "\t" << colBoard << std::endl;
0215         mapping[chCAEN][0]    = layerMod; 
0216         mapping[chCAEN][1]    = chBoard; 
0217         mapping[chCAEN][2]    = rowBoard; 
0218         mapping[chCAEN][3]    = colBoard; 
0219         nChmapped++;
0220         backwardMapping[chBoard][layerMod] = chCAEN;
0221         backwardMappingBoard[chBoard][layerMod] = ((TString)((TObjString*)tempArr->At(0))->GetString()).Atoi();
0222         if (verbosity > 0) std::cout << backwardMapping[chBoard][layerMod] << std::endl;
0223         delete tempArr;
0224     }
0225     
0226     for (Int_t ch = 0; ch< 64; ch++){
0227       if (verbosity > 0) std::cout  << "channel: " << ch << " location plane: " << mapping[ch][0] << "\t tile: " << mapping[ch][1] << std::endl;
0228     }
0229   
0230   
0231     Int_t nActiveLayers   = 0;
0232     Int_t minActiveLayer  = -1;
0233     Int_t maxActiveLayer  = -1;
0234     Int_t nActiveRBCh     = 0;
0235     Int_t minActiveRBCh   = -1;
0236     Int_t maxActiveRBCh   = -1;
0237 
0238     std::cout << "===================================================================" << std::endl;
0239     std::cout << "================ Visualization of layers ==========================" << std::endl;
0240     std::cout << "===================================================================" << std::endl;
0241     for (Int_t l = 0; l < gMaxLayers; l++){
0242         for (Int_t c = 1; c < 9; c++){
0243             std::cout << backwardMapping[c][l] << "\t || \t" ;
0244             if (lActive[l] == kFALSE ){
0245                 if ( backwardMapping[c][l] != -1 ) lActive[l] = kTRUE;
0246             }
0247             if (lActiveCh[c] == kFALSE ){
0248               if ( backwardMapping[c][l] != -1 ) lActiveCh[c] = kTRUE;
0249             }
0250         }
0251         std::cout << std::endl;
0252     }
0253     std::cout << "===================================================================" << std::endl;
0254     
0255     // count total number of active layers according to mapping file and determine maximum layer
0256     for (Int_t l = 0; l<gMaxLayers; l++){
0257       if (lActive[l]){
0258         nActiveLayers++; 
0259         maxActiveLayer = l;
0260         if (minActiveLayer == -1) minActiveLayer = l;
0261       }
0262     }
0263     // count total number of active channel according to mapping file and determine maximum RBchannel
0264     for (Int_t c = 0; c<9; c++){
0265       if (lActiveCh[c]){
0266         nActiveRBCh++; 
0267         maxActiveRBCh = c;
0268         if (minActiveRBCh == -1) minActiveRBCh = c;
0269       }
0270     }
0271     std::cout << "There are " << nActiveLayers << " active layers. The lowest layer is " << minActiveLayer << "  and the highest layer is " << maxActiveLayer << std::endl;
0272     std::cout << "There are " << nActiveRBCh << " active read-outboard channels. The lowest channel is " << minActiveRBCh << "  and the highest channel is " << maxActiveRBCh << std::endl;
0273     
0274     // ********************************************************************************************************
0275     // RAW data monitoring histos per CAEN board and channel
0276     // ********************************************************************************************************    
0277     TH1D* histLG[gMaxBoard][gMaxChannels];                      // low gain all triggers
0278     TH1D* histHG[gMaxBoard][gMaxChannels];                      // high gain all triggers
0279     TH2D* hist_T_HG[gMaxBoard][gMaxChannels];                   // time dependent high gain all triggers
0280     TH2D* histLGHG[gMaxBoard][gMaxChannels];                    // LG - HG correlation
0281     TH2D* histHGLG[gMaxBoard][gMaxChannels];                    // HG - LG correlation
0282     //***************************************************
0283     // mapped channels - layer & read-out board channels
0284     //***************************************************
0285     // channel counting in layer starts with 1, easier acces make array go to 8
0286     TH1D* histLG_mapped[gMaxLayers][9];                         // low gain all triggers
0287     TH1D* histHG_mapped[gMaxLayers][9];                         // high gain all triggers
0288     
0289     //***************************************************
0290     // map of triggers above threshold
0291     //***************************************************
0292     TH3D* hist3DMap = new TH3D("h_map_z_x_y", "; layer; col; row", 14, -0.5, 13.5, 4, -0.5, 3.5, 2, -0.5, 1.5);
0293     TH2D* hist2DMap = new TH2D("h_map_z_channel", "; channel; layer", 8, 0.5, 8.5, 14, -0.5, 13.5);
0294     TH1D* hist1DMap = new TH1D("h_map_channel", "; channel", 8, 0.5, 8.5);
0295     
0296     //***************************************************
0297     TH1D* histNChAboveNoise   = new TH1D("h_NchannelAboveNoise", "; N_{channels}; counts", (gMaxBoard*gMaxChannels+1), -0.5, (gMaxBoard*gMaxChannels-0.5));
0298     
0299     //***************************************************
0300     // Init histograms
0301     //***************************************************
0302     for (Int_t j = 0; j < gMaxBoard; j++){
0303       for (Int_t i = 0; i < gMaxChannels; i++){
0304         hist_T_HG[j][i]    = new TH2D(Form("h_T_HG_B%d_C%02d",j,i),"; t (min); HG (adc); counts",2000,0,100,4002,-1,2000);
0305         histHG[j][i]    = new TH1D(Form("h_HG_B%d_C%02d",j,i),"; HG (adc); counts",4201,-200,4001);
0306         histLG[j][i]    = new TH1D(Form("h_LG_B%d_C%02d",j,i),"; LG (adc); counts",4201,-200,4001);
0307         histLGHG[j][i]  = new TH2D(Form("h_LGHG_B%d_C%02d",j,i),"; LG (adc); HG (adc)",4200/5,-200,4001,4200/5,-200,4001);
0308         histHGLG[j][i]  = new TH2D(Form("h_HGLG_B%d_C%02d",j,i),"; LG (adc); HG (adc)",4200/5,-200,4001,4200/5,-200,4001);
0309       }
0310     }
0311     for (Int_t l = 0; l < gMaxLayers; l++){
0312       for (Int_t c = 0; c < 9; c++){
0313         histLG_mapped[l][c]             = nullptr;
0314         histHG_mapped[l][c]             = nullptr; 
0315       }
0316     }
0317     
0318     // ********************************************************************************************************
0319     // Setup of global variables for parsing tree
0320     // ********************************************************************************************************
0321     Long64_t adcThreshold = 100;                    // minimum ADC count for triggers
0322     Double_t scaledThr        = adcThreshold*1.5;   // scaled ADC trigger threshold
0323     Double_t tstapMin         = 0;                  // min time stamp after start of trigger (#mus)
0324     Double_t tstapMax         = 0;                  // current max time stamp after start of trigger (#mus)
0325     
0326     Long64_t nEventsProcessed = 0;              // running counter
0327     // use this if you wanna start at a specific event for debug
0328     Long64_t startEvent = 0;
0329     if (minNEvent > 0) startEvent = minNEvent;
0330     
0331     // ********************************************************************************************************
0332     // First loop over full tree to obtain RAW histograms
0333     // ********************************************************************************************************
0334     for (Long64_t i=startEvent; i<nEntriesTree;i++) {
0335         // load current event
0336         tt_event->GetEntry(i);
0337         if (i == startEvent) tstapMin = gTRtimeStamp;
0338         tstapMax = gTRtimeStamp;
0339         Double_t tCurr= (tstapMax-tstapMin)/1e6/60; // elapsed time in min
0340         nEventsProcessed++;
0341 
0342         Int_t nChNoNoise = 0;             // number of channels above noise level
0343         // processing progress info
0344         if(i>0 && nEntriesTree>100 && i%(nEntriesTree/(20))==0) std::cout << "//processed " << 100*(i)/nEntriesTree << "%"  << std::endl;
0345         if(verbosity>1){
0346           std::cout << "***********************************************************************************************************" << std::endl;
0347           std::cout << "event " << i << std::endl;
0348           std::cout << "***********************************************************************************************************" << std::endl;
0349         }
0350         // temporary mapping of channels
0351         Int_t signal[2][64][3] = {{{0}}}; // [board][channel][HG, LG, trigg HG]
0352         // readEvent
0353         if (verbosity > 1)std::cout << "------------- Event -------------------------" << std::endl;
0354         if (verbosity > 1)std::cout << gTrID << "\t" << gTRtimeStamp << std::endl;
0355         for(Int_t ch=0; ch<gMaxChannels && ch < nChmapped ; ch++){   
0356           signal[gBoard[ch]][gChannel[ch]][0] = gHG[ch];
0357           signal[gBoard[ch]][gChannel[ch]][1] = gLG[ch];
0358           if (gHG[ch] > scaledThr)  // fill potential trigger signals
0359             signal[gBoard[ch]][gChannel[ch]][2] = gHG[ch]-scaledThr;
0360           
0361           // filling hists
0362           histHG[gBoard[ch]][gChannel[ch]]->Fill(gHG[ch]); 
0363           hist_T_HG[gBoard[ch]][gChannel[ch]]->Fill(tCurr,gHG[ch]); 
0364           histLG[gBoard[ch]][gChannel[ch]]->Fill(gLG[ch]); 
0365           histLGHG[gBoard[ch]][gChannel[ch]]->Fill(gLG[ch], gHG[ch]); 
0366           histHGLG[gBoard[ch]][gChannel[ch]]->Fill(gHG[ch], gLG[ch]); 
0367           if (verbosity > 1)std::cout << "--> board: "<< gBoard[ch] << "\t ch:"<< gChannel[ch] << "\t LG:" << gLG[ch] << "\t HG:" << gHG[ch] << std::endl;
0368           if (gHG[ch] > adcThreshold){
0369             nChNoNoise++;
0370             if (verbosity > 1) std::cout << "not noise" << std::endl;
0371             hist3DMap->Fill(mapping[gChannel[ch]][0],mapping[gChannel[ch]][3],mapping[gChannel[ch]][2]);
0372             hist2DMap->Fill(mapping[gChannel[ch]][1],mapping[gChannel[ch]][0]);
0373             hist1DMap->Fill(mapping[gChannel[ch]][1]);
0374           }
0375         }
0376         // fill channels above noise
0377         histNChAboveNoise->Fill(nChNoNoise);
0378         if (verbosity > 1) std::cout << "Channels above noise: "<< nChNoNoise << std::endl;        
0379     }
0380     
0381     //**********************************************************************
0382     // Monitoring infos
0383     //**********************************************************************
0384     std::cout << "Total events processed:" << nEventsProcessed << std::endl;
0385     Double_t tdiff = (tstapMax-tstapMin)/1e6/60; //time in min
0386     std::cout << "times: " << tstapMin << "\t" << tstapMax <<"\t elapsed: " <<  tdiff  << " min"<< std::endl;
0387     
0388     //**********************************************************************
0389     // Create canvases for single channel plotting
0390     //**********************************************************************
0391     TCanvas* canvas2DCorr = new TCanvas("canvasCorrPlots","",0,0,1350,1200);  // gives the page size
0392     DefaultCancasSettings( canvas2DCorr, 0.1, 0.1, 0.02, 0.08);
0393     canvas2DCorr->SetLogz();
0394     TCanvas* canvas1DNoise = new TCanvas("canvas1DNoise","",0,0,700,500);  // gives the page size
0395     DefaultCancasSettings( canvas1DNoise, 0.07, 0.02, 0.02, 0.08);
0396     canvas1DNoise->SetLogy();
0397     TCanvas* canvas1DDiffTrigg = new TCanvas("canvas1DDiffTrigg","",0,0,700,500);  // gives the page size
0398     DefaultCancasSettings( canvas1DDiffTrigg, 0.07, 0.02, 0.02, 0.08);
0399     canvas1DDiffTrigg->SetLogy();
0400     
0401     //***********************************************************************************************************
0402     //********************************* 8 Panel overview plot  **************************************************
0403     //***********************************************************************************************************
0404     //*****************************************************************
0405     // Test beam geometry (beam coming from viewer)
0406     //==============================================
0407     //||    8    ||    7    ||    6    ||    5    ||
0408     //==============================================
0409     //||    1    ||    2    ||    3    ||    4    ||
0410     //==============================================
0411     // rebuild pad geom in similar way (numbering -1)
0412     //*****************************************************************
0413     TCanvas* canvas8Panel;
0414     TPad* pad8Panel[8];
0415     Double_t topRCornerX[8];
0416     Double_t topRCornerY[8];
0417     Int_t textSizePixel = 30;
0418     Double_t relSize8P[8];
0419     CreateCanvasAndPadsFor8PannelTBPlot(canvas8Panel, pad8Panel,  topRCornerX, topRCornerY, relSize8P, textSizePixel);
0420     
0421     //**********************************************************************
0422     // Fits for noise and LG/HG calibration
0423     //**********************************************************************
0424     TF1* fitGausHG_BG[gMaxBoard][gMaxChannels];
0425     TF1* fitGausLG_BG[gMaxBoard][gMaxChannels];
0426     TF1* fitGausHG_BG_mapped[gMaxLayers][9];
0427     TF1* fitGausLG_BG_mapped[gMaxLayers][9];
0428     TF1* fitLGHGCorr[gMaxBoard][gMaxChannels];
0429     TF1* fitHGLGCorr[gMaxBoard][gMaxChannels];
0430     Double_t mean[4][gMaxBoard][gMaxChannels];
0431     Double_t sigma[4][gMaxBoard][gMaxChannels];
0432     Double_t cslope[4][gMaxBoard][gMaxChannels];
0433     Double_t coffset[4][gMaxBoard][gMaxChannels];
0434 
0435     //**********************************************************************
0436     // Monitoring hists for fits to noise and slope LG-HG
0437     //**********************************************************************
0438     TH1D* histNoiseSigma_HG     = new TH1D("histNoiseSigma_HG", "; noise #sigma (HG ADC)", 400, 0, 20);
0439     TH1D* histNoiseMean_HG      = new TH1D("histNoiseMean_HG", "; noise #mu (HG ADC)", 400, 30, 70);
0440     TH1D* histNoiseSigma_LG     = new TH1D("histNoiseSigma_LG", "; noise #sigma (LG ADC)", 400, 0, 20);
0441     TH1D* histNoiseMean_LG      = new TH1D("histNoiseMean_LG", "; noise #mu (LG ADC)", 400, 30, 70);
0442     TH1D* histLGHG_slope        = new TH1D("histLGHGslope", "; slope (HG adc/LG adc)", 400, 0, 40);
0443     
0444     TH2D* hist2DNoiseSigma_HG   = new TH2D("hist2DNoiseSigma_HG_z_channel", "; channel; layer; noise #sigma (HG ADC)", 8, 0.5, 8.5, 14, -0.5, 13.5);
0445     TH2D* hist2DNoiseMean_HG    = new TH2D("hist2DNoiseMean_HG_z_channel", "; channel; layer; noise #mu (HG ADC)", 8, 0.5, 8.5, 14, -0.5, 13.5);
0446     TH2D* hist2DNoiseSigma_LG   = new TH2D("hist2DNoiseSigma_LG_z_channel", "; channel; layer; noise #sigma (LG ADC)", 8, 0.5, 8.5, 14, -0.5, 13.5);
0447     TH2D* hist2DNoiseMean_LG    = new TH2D("hist2DNoiseMean_LG_z_channel", "; channel; layer; noise #mu (LG ADC)", 8, 0.5, 8.5, 14, -0.5, 13.5);
0448     TH2D* hist2DLGHG_slope      = new TH2D("hist2D_LGHGslope_z_channel", "; channel; layer; slope (HG adc/LG adc)", 8, 0.5, 8.5, 14, -0.5, 13.5);
0449     TH2D* hist2DLGHG_offset     = new TH2D("hist2D_LGHGoffset_z_channel", "; channel; layer; offset (HG ADC)", 8, 0.5, 8.5, 14, -0.5, 13.5);
0450     
0451     // 1D channel representation of fit values, x axis scales as 10x layer count + channel within one assembley, 
0452     // - layer 3 channel 3: 33
0453     // - layer 0 channel 2: 2
0454     Int_t tempLayerBinning = maxActiveLayer;
0455     if (maxActiveLayer == 0) tempLayerBinning = 1;
0456     TH1D* hist1DNoiseSigma_HG   = new TH1D("hist1DNoiseSigma_HG_channels", "; 10x layer + board channel; noise #sigma (HG ADC)", 10*tempLayerBinning+1, -0.5, 10*tempLayerBinning+0.5 );
0457     TH1D* hist1DNoiseMean_HG    = new TH1D("hist1DNoiseMean_HG_channels", "; 10x layer + board channel; noise #mu (HG ADC)", 10*tempLayerBinning+1, -0.5, 10*tempLayerBinning+0.5 );
0458     TH1D* hist1DNoiseSigma_LG   = new TH1D("hist1DNoiseSigma_LG_channels", "; 10x layer + board channel; noise #sigma (LG ADC)", 10*tempLayerBinning+1, -0.5, 10*tempLayerBinning+0.5 );
0459     TH1D* hist1DNoiseMean_LG    = new TH1D("hist1DNoiseMean_LG_channels", "; 10x layer + board channel; noise #mu (LG ADC)", 10*tempLayerBinning+1, -0.5, 10*tempLayerBinning+0.5 );
0460     TH1D* hist1DLGHG_slope      = new TH1D("hist1DLGHG_slope_channels", "; 10x layer + board channel; slope (HG adc/LG adc)", 10*tempLayerBinning+1, -0.5, 10*tempLayerBinning+0.5 );
0461     TH1D* hist1DLGHG_offset     = new TH1D("hist1DLGHG_slope_channels", "; 10x layer + board channel; offset (HG adc/LG adc)", 10*tempLayerBinning+1, -0.5, 10*tempLayerBinning+0.5 );
0462     TH1D* hist1DHGLG_slope      = new TH1D("hist1DHGLG_slope_channels", "; 10x layer + board channel; slope (LG adc/HG adc)", 10*tempLayerBinning+1, -0.5, 10*tempLayerBinning+0.5 );
0463     TH1D* hist1DHGLG_offset     = new TH1D("hist1DHGLG_slope_channels", "; 10x layer + board channel; offset (LG adc/HG adc)", 10*tempLayerBinning+1, -0.5, 10*tempLayerBinning+0.5 );
0464     
0465     // 1D channel representation of fit values, x axis scales as CAEN board channelns 64x CAEN board # + CAEN channel
0466     TH1D* hist1DCAEN_NoiseSigma_HG   = new TH1D("hist1DCAEN_NoiseSigma_HG_channels", "; 64x CAEN board + CAEN channel; noise #sigma (HG ADC)", 64*gMaxBoard+1, -0.5, 64*gMaxBoard+0.5 );
0467     TH1D* hist1DCAEN_NoiseMean_HG    = new TH1D("hist1DCAEN_NoiseMean_HG_channels", "; 64x CAEN board + CAEN channel; noise #mu (HG ADC)", 64*gMaxBoard+1, -0.5, 64*gMaxBoard+0.5 );
0468     TH1D* hist1DCAEN_NoiseSigma_LG   = new TH1D("hist1DCAEN_NoiseSigma_LG_channels", "; 64x CAEN board + CAEN channel; noise #sigma (LG ADC)", 64*gMaxBoard+1, -0.5, 64*gMaxBoard+0.5 );
0469     TH1D* hist1DCAEN_NoiseMean_LG    = new TH1D("hist1DCAEN_NoiseMean_LG_channels", "; 64x CAEN board + CAEN channel; noise #mu (LG ADC)", 64*gMaxBoard+1, -0.5, 64*gMaxBoard+0.5 );
0470     TH1D* hist1DCAEN_LGHG_slope      = new TH1D("hist1DCAEN_LGHG_slope_channels", "; 64x CAEN board + CAEN channel; slope (HG adc/LG adc)", 64*gMaxBoard+1, -0.5, 64*gMaxBoard+0.5 );
0471     TH1D* hist1DCAEN_LGHG_offset     = new TH1D("hist1DCAEN_LGHG_slope_channels", "; 64x CAEN board + CAEN channel; offset (HG adc/LG adc)", 64*gMaxBoard+1, -0.5, 64*gMaxBoard+0.5 );
0472     TH1D* hist1DCAEN_HGLG_slope      = new TH1D("hist1DCAEN_HGLG_slope_channels", "; 64x CAEN board + CAEN channel; slope (LG adc/HG adc)", 64*gMaxBoard+1, -0.5, 64*gMaxBoard+0.5 );
0473     TH1D* hist1DCAEN_HGLG_offset     = new TH1D("hist1DCAEN_HGLG_slope_channels", "; 64x CAEN board + CAEN channel; offset (LG adc/HG adc)", 64*gMaxBoard+1, -0.5, 64*gMaxBoard+0.5 );
0474 
0475     //**********************************************************************
0476     // Initialize fits for all layers and channels to nullptr
0477     //**********************************************************************    
0478     for (Int_t l = 0; l < gMaxLayers; l++){
0479       for (Int_t c = 0; c< 9;c++){
0480         fitGausHG_BG_mapped[l][c] = nullptr;
0481         fitGausLG_BG_mapped[l][c] = nullptr;
0482       }
0483     }
0484     
0485     // ********************************************************************************************************
0486     // Fitting Noise and correlation HG-LG, LG-HG
0487     // -> detailed plotting per channel
0488     // -> remapping of CAEN board+channels to physicsal layers and readout-channels
0489     // ********************************************************************************************************    
0490     for (Int_t j = 0; j < gMaxBoard; j++){
0491       for (Int_t i = 0; i < gMaxChannels; i++){
0492         mean[0][j][i] = -1;
0493         mean[1][j][i] = -1;
0494         mean[2][j][i] = -1;
0495         mean[3][j][i] = -1;
0496         sigma[0][j][i] = -1;
0497         sigma[1][j][i] = -1;
0498         sigma[2][j][i] = -1;
0499         sigma[3][j][i] = -1;
0500 
0501         if (! ((histHG[j][i] && histHG[j][i]->GetEntries() > 0)|| (histLG[j][i] && histLG[j][i]->GetEntries() > 0)) ){
0502           fitGausLG_BG[j][i]  = nullptr;
0503           fitGausHG_BG[j][i]  = nullptr;
0504           fitLGHGCorr[j][i]   = nullptr;
0505           fitHGLGCorr[j][i]   = nullptr;
0506           continue;
0507         }
0508         Int_t chMap     = j*64 + i;
0509         Int_t chBoard   = mapping[chMap][1];
0510         Int_t layer     = mapping[chMap][0];
0511         Int_t channelBin1D = hist1DNoiseSigma_HG->FindBin(layer*10+chBoard);          
0512           
0513         if (verbosity > 0)std::cout << j << "\t" << i << "\t" << chMap << "\t L: " << layer  << "\t C:" <<  chBoard << "\t bin ID "<< channelBin1D << std::endl;
0514         
0515         // ********************************************************
0516         // map raw and trigger histos from CAEN numbering to 
0517         // physical layers & channels
0518         // ********************************************************
0519         histHG_mapped[layer][chBoard] = (TH1D*)histHG[j][i]->Clone(Form("h_HG_mapped_L%d_C%02d",layer,chBoard));
0520         histLG_mapped[layer][chBoard] = (TH1D*)histLG[j][i]->Clone(Form("h_LG_mapped_L%d_C%02d",layer,chBoard));
0521         // ********************************************************
0522         // Fit & plot noise histos
0523         // -> map noise histos from CAEN numb. -> physical
0524         // -> fill monitoring plots
0525         // ********************************************************
0526         Bool_t bFit = kFALSE;
0527         // HG noise fits
0528         bFit = FitNoiseWithBG (histHG[j][i], fitGausHG_BG[j][i], mean[0][j][i], mean[1][j][i], sigma[0][j][i], sigma[1][j][i], j, i, "f_GaussBG_HG", "HG", verbosity);
0529         if (bFit && bDetPlot) PlotNoiseSingle (canvas1DNoise, histHG[j][i], fitGausHG_BG[j][i], mean[0][j][i], mean[1][j][i], sigma[0][j][i], sigma[1][j][i], j, i, layer, chBoard, 
0530                                                 Form("%s/HG_NoiseWithFit", outputDirPlotsDet.Data()), currentRunInfo, 0.04);
0531         if (bFit){
0532           fitGausHG_BG_mapped[layer][chBoard] = (TF1*)fitGausHG_BG[j][i]->Clone(Form("f_GaussBG_HG_mapped_L%d_C%02d",layer,chBoard));
0533         }
0534         if (bFit == kFALSE) fitGausHG_BG[j][i] = nullptr;
0535         // fill fit monitoring hists HG
0536         if (bFit){
0537           hist2DNoiseMean_HG->Fill(chBoard,layer,mean[0][j][i]);
0538           hist2DNoiseSigma_HG->Fill(chBoard,layer,sigma[0][j][i]);
0539           histNoiseMean_HG->Fill(mean[0][j][i]);
0540           histNoiseSigma_HG->Fill(sigma[0][j][i]);
0541           
0542           hist1DNoiseSigma_HG->SetBinContent(channelBin1D, sigma[0][j][i]);
0543           hist1DNoiseSigma_HG->SetBinError(channelBin1D, sigma[1][j][i]);
0544           hist1DNoiseMean_HG->SetBinContent(channelBin1D, mean[0][j][i]);
0545           hist1DNoiseMean_HG->SetBinError(channelBin1D, mean[1][j][i]);
0546           hist1DCAEN_NoiseSigma_HG->SetBinContent(chMap, sigma[0][j][i]);
0547           hist1DCAEN_NoiseSigma_HG->SetBinError(chMap, sigma[1][j][i]);
0548           hist1DCAEN_NoiseMean_HG->SetBinContent(chMap, mean[0][j][i]);          
0549           hist1DCAEN_NoiseMean_HG->SetBinError(chMap, mean[1][j][i]);          
0550         }
0551         // reset boolean for fit success monitoring
0552         bFit = kFALSE;
0553         // LG noise fits
0554         bFit = FitNoiseWithBG (histLG[j][i], fitGausLG_BG[j][i], mean[2][j][i], mean[3][j][i], sigma[2][j][i], sigma[3][j][i], j, i, "f_GaussBG_LG", "LG");
0555         if (bFit && bDetPlot) PlotNoiseSingle (canvas1DNoise, histLG[j][i], fitGausLG_BG[j][i], mean[2][j][i], mean[3][j][i], sigma[2][j][i], sigma[3][j][i], j, i, layer, chBoard, 
0556                                                 Form("%s/LG_NoiseWithFit", outputDirPlotsDet.Data()), currentRunInfo, 0.04);
0557         if (bFit){
0558           fitGausLG_BG_mapped[layer][chBoard] = (TF1*)fitGausLG_BG[j][i]->Clone(Form("f_GaussBG_LG_mapped_L%d_C%02d",layer,chBoard));
0559         }
0560         if (bFit == kFALSE) fitGausLG_BG[j][i] = nullptr;
0561         // fill fit monitoring hists LG
0562         if (bFit){
0563           hist2DNoiseMean_LG->Fill(chBoard,layer,mean[2][j][i]);
0564           hist2DNoiseSigma_LG->Fill(chBoard,layer,sigma[2][j][i]);
0565           histNoiseMean_LG->Fill(mean[2][j][i]);
0566           histNoiseSigma_LG->Fill(sigma[2][j][i]);
0567           hist1DNoiseSigma_LG->SetBinContent(channelBin1D, sigma[2][j][i]);
0568           hist1DNoiseSigma_LG->SetBinError(channelBin1D, sigma[3][j][i]);
0569           hist1DNoiseMean_LG->SetBinContent(channelBin1D, mean[2][j][i]);
0570           hist1DNoiseMean_LG->SetBinError(channelBin1D, mean[3][j][i]);
0571           hist1DCAEN_NoiseSigma_LG->SetBinContent(chMap, sigma[2][j][i]);
0572           hist1DCAEN_NoiseSigma_LG->SetBinError(chMap, sigma[3][j][i]);
0573           hist1DCAEN_NoiseMean_LG->SetBinContent(chMap, mean[2][j][i]);          
0574           hist1DCAEN_NoiseMean_LG->SetBinError(chMap, mean[3][j][i]);          
0575         }
0576         // ********************************************************
0577         // Plot 2D distribution time vs HG for monitoring
0578         // ********************************************************
0579         if (hist_T_HG[j][i]  && hist_T_HG[j][i]->GetEntries() > 0 && bDetPlot){
0580           canvas2DCorr->cd();
0581             SetStyleHistoTH2ForGraphs( hist_T_HG[j][i], hist_T_HG[j][i]->GetXaxis()->GetTitle(), hist_T_HG[j][i]->GetYaxis()->GetTitle(), 0.85*textSizeRel, textSizeRel, 0.85*textSizeRel, textSizeRel,0.9, 1.25);  
0582             // find max x bin
0583             hist_T_HG[j][i]->GetXaxis()->SetRangeUser(0,tdiff+0.1);
0584             hist_T_HG[j][i]->Draw("colz");
0585             
0586             TLatex *labelChannel                     = new TLatex(0.85,0.92,Form("CAEN: B %d, C %d, Stack: L %d, C%d",j, i, layer, chBoard));
0587             SetStyleTLatex( labelChannel, textSizeRel,4,1,42,kTRUE,31);
0588             labelChannel->Draw();
0589 
0590           canvas2DCorr->SaveAs(Form("%s/T_HG_B%d_C%02d.pdf", outputDirPlotsDet.Data(), j,i));
0591           delete labelChannel;
0592         }
0593         // ********************************************************
0594         // Fit & Plot 2D LG vs HG
0595         // ********************************************************
0596         if (histLGHG[j][i] && histLGHG[j][i]->GetEntries() > 0){
0597           FitAndPlotGainCorr ( histLGHG[j][i], fitLGHGCorr[j][i], "f_LGHGCorr", 100, 380, 500, 4000,
0598                               cslope[0][j][i], cslope[1][j][i], coffset[0][j][i], coffset[1][j][i], 
0599                               j, i, layer, chBoard,
0600                               kTRUE, canvas2DCorr, Form("%s/LG_HG_Corr", outputDirPlotsDet.Data()), textSizeRel);
0601           hist2DLGHG_slope->Fill(chBoard,layer,cslope[0][j][i]);
0602           hist2DLGHG_offset->Fill(chBoard,layer,coffset[0][j][i]);
0603           histLGHG_slope->Fill(cslope[0][j][i]);
0604           hist1DLGHG_slope->SetBinContent(channelBin1D, cslope[0][j][i]);
0605           hist1DLGHG_slope->SetBinError(channelBin1D, cslope[1][j][i]);
0606           hist1DLGHG_offset->SetBinContent(channelBin1D, coffset[0][j][i]);
0607           hist1DLGHG_offset->SetBinError(channelBin1D, coffset[1][j][i]);
0608           hist1DCAEN_LGHG_slope->SetBinContent(chMap, cslope[0][j][i]);
0609           hist1DCAEN_LGHG_slope->SetBinError(chMap, cslope[1][j][i]);
0610           hist1DCAEN_LGHG_offset->SetBinContent(chMap, coffset[0][j][i]);
0611           hist1DCAEN_LGHG_offset->SetBinError(chMap, coffset[1][j][i]);
0612         } else {
0613           fitLGHGCorr[j][i] = nullptr;
0614         }
0615         if (histHGLG[j][i] && histLGHG[j][i]->GetEntries() > 0){
0616           FitAndPlotGainCorr ( histHGLG[j][i], fitHGLGCorr[j][i], "f_HGLGCorr", 100, 3800, 4000, 500,
0617                               cslope[2][j][i], cslope[3][j][i], coffset[2][j][i], coffset[3][j][i], 
0618                               j, i, layer, chBoard,
0619                               kFALSE, canvas2DCorr, Form("%s/HG_LG_Corr", outputDirPlotsDet.Data()), textSizeRel);
0620           hist1DHGLG_slope->SetBinContent(channelBin1D, cslope[2][j][i]);
0621           hist1DHGLG_slope->Fill(channelBin1D, cslope[3][j][i]);
0622           hist1DHGLG_offset->SetBinContent(channelBin1D, coffset[2][j][i]);
0623           hist1DHGLG_offset->SetBinError(channelBin1D, coffset[3][j][i]);
0624           hist1DCAEN_HGLG_slope->SetBinContent(chMap, cslope[2][j][i]);
0625           hist1DCAEN_HGLG_slope->SetBinError(chMap, cslope[3][j][i]);
0626           hist1DCAEN_HGLG_offset->SetBinContent(chMap, coffset[2][j][i]);
0627           hist1DCAEN_HGLG_offset->SetBinError(chMap, coffset[3][j][i]);
0628         } else {
0629           fitHGLGCorr[j][i] = nullptr;
0630         }        
0631         
0632         // ********************************************************
0633         // Plot different triggers together
0634         // ********************************************************
0635         Int_t hgmax = 1050;
0636         std::cout << __LINE__ << std::endl;
0637         if (bDetPlot)PlotOverlayDiffTriggers( canvas1DDiffTrigg, histHG[j][i], nullptr, nullptr, fitGausHG_BG[j][i],
0638                                               0, hgmax, Form("%s/HG_DiffTriggers", outputDirPlotsDet.Data()),
0639                                               j, i, layer, chBoard, currentRunInfo, 0.04);
0640         std::cout << __LINE__ << std::endl;
0641         if (bDetPlot)PlotOverlayDiffTriggers( canvas1DDiffTrigg, histLG[j][i], nullptr, nullptr, fitGausLG_BG[j][i],
0642                                               0, 1050, Form("%s/LG_DiffTriggers", outputDirPlotsDet.Data()),
0643                                               j, i, layer, chBoard, currentRunInfo, 0.04);
0644       }
0645     }
0646     
0647     for (Int_t l = 0; l < gMaxLayers; l++){
0648       if (!lActive[l]) continue;     
0649       std::cout << __LINE__ << std::endl;
0650       PlotDiffTriggersFullLayer (canvas8Panel,pad8Panel, topRCornerX, topRCornerY, relSize8P, textSizePixel, 
0651                                 histHG_mapped[l], nullptr, nullptr, fitGausHG_BG_mapped[l], 
0652                                 0, 1000, 1.2, l , Form("%s/TriggerOverlay_HG_Zoomed_Layer%02d.pdf" ,outputDirPlots.Data(), l), currentRunInfo);
0653       std::cout << __LINE__ << std::endl;
0654       PlotDiffTriggersFullLayer (canvas8Panel,pad8Panel, topRCornerX, topRCornerY, relSize8P, textSizePixel, 
0655                                 histLG_mapped[l], nullptr, nullptr, fitGausLG_BG_mapped[l], 
0656                                 0, 800, 1.2, l , Form("%s/TriggerOverlay_LG_Zoomed_Layer%02d.pdf" ,outputDirPlots.Data(), l), currentRunInfo);
0657     // ********************************************************************************************************
0658     // Overlay in same layer
0659     // ********************************************************************************************************
0660       std::cout << __LINE__ << std::endl;
0661       PlotChannelOverlaySameLayer( canvas1DDiffTrigg, histHG_mapped[l], 1, 9,
0662                                   -100, 1000, 1./5, Form("%s/HG", outputDirPlots.Data()), l, currentRunInfo, 0.04,"hist");
0663       std::cout << __LINE__ << std::endl;
0664       PlotChannelOverlaySameLayerWithFitsBG( canvas1DDiffTrigg, histHG_mapped[l], fitGausHG_BG_mapped[l], 1, 9,
0665                                   -100, 1000, 1./5, Form("%s/HGWithFits", outputDirPlots.Data()), l, currentRunInfo, 0.04,"hist");
0666       std::cout << __LINE__ << std::endl;
0667       PlotChannelOverlaySameLayer( canvas1DDiffTrigg, histLG_mapped[l], 1, 9,
0668                                   -100, 500, 1, Form("%s/LG", outputDirPlots.Data()), l, currentRunInfo, 0.04);        
0669     
0670     }    
0671     
0672     // ********************************************************
0673     // Create graphs per board and channels with calib values
0674     // -> noise mean & sigma LG, HG
0675     // -> slope and offset of LG-HG & HG-LG correlation
0676     // ********************************************************
0677     TGraphErrors* gNoiseMeanHG      = CreateGraphFromHistAndCleanup(hist1DNoiseMean_HG, "graph_mean_Noise_HG_channels");
0678     TGraphErrors* gNoiseMeanLG      = CreateGraphFromHistAndCleanup(hist1DNoiseMean_LG, "graph_mean_Noise_LG_channels");
0679     TGraphErrors* gNoiseSigmaHG     = CreateGraphFromHistAndCleanup(hist1DNoiseSigma_HG, "graph_sigma_Noise_HG_channels");
0680     TGraphErrors* gNoiseSigmaLG     = CreateGraphFromHistAndCleanup(hist1DNoiseSigma_LG, "graph_sigma_Noise_LG_channels");
0681     TGraphErrors* gCorrLGHGSlope    = CreateGraphFromHistAndCleanup(hist1DLGHG_slope, "graph_slope_corr_LGHG_channels");
0682     TGraphErrors* gCorrLGHGOffset   = CreateGraphFromHistAndCleanup(hist1DLGHG_offset, "graph_offset_corr_LGHG_channels");
0683     TGraphErrors* gCorrHGLGSlope    = CreateGraphFromHistAndCleanup(hist1DHGLG_slope, "graph_slope_corr_HGLG_channels");
0684     TGraphErrors* gCorrHGLGOffset   = CreateGraphFromHistAndCleanup(hist1DHGLG_offset, "graph_offset_corr_HGLG_channels");
0685 
0686     TGraphErrors* gCAEN_NoiseMeanHG     = CreateGraphFromHistAndCleanup(hist1DCAEN_NoiseMean_HG, "graphCAEN_mean_Noise_HG_channels");
0687     TGraphErrors* gCAEN_NoiseMeanLG     = CreateGraphFromHistAndCleanup(hist1DCAEN_NoiseMean_LG, "graphCAEN_mean_Noise_LG_channels");
0688     TGraphErrors* gCAEN_NoiseSigmaHG    = CreateGraphFromHistAndCleanup(hist1DCAEN_NoiseSigma_HG, "graphCAEN_sigma_Noise_HG_channels");
0689     TGraphErrors* gCAEN_NoiseSigmaLG    = CreateGraphFromHistAndCleanup(hist1DCAEN_NoiseSigma_LG, "graphCAEN_sigma_Noise_LG_channels");
0690     TGraphErrors* gCAEN_CorrLGHGSlope   = CreateGraphFromHistAndCleanup(hist1DCAEN_LGHG_slope, "graphCAEN_slope_corr_LGHG_channels");
0691     TGraphErrors* gCAEN_CorrLGHGOffset  = CreateGraphFromHistAndCleanup(hist1DCAEN_LGHG_offset, "graphCAEN_offset_corr_LGHG_channels");
0692     TGraphErrors* gCAEN_CorrHGLGSlope   = CreateGraphFromHistAndCleanup(hist1DCAEN_HGLG_slope, "graphCAEN_slope_corr_HGLG_channels");
0693     TGraphErrors* gCAEN_CorrHGLGOffset  = CreateGraphFromHistAndCleanup(hist1DCAEN_HGLG_offset, "graphCAEN_offset_corr_HGLG_channels");
0694         
0695     // ********************************************************************************************************
0696     // Noise subtracted data monitoring histos per CAEN board and channel
0697     // ********************************************************************************************************    
0698     TH1D* histNSLG[gMaxBoard][gMaxChannels];                // low gain all triggers
0699     TH1D* histNSHG[gMaxBoard][gMaxChannels];                // high gain all triggers
0700     TH1D* histNSHG_BG[gMaxBoard][gMaxChannels];             // high gain all triggers backgrounds
0701     TH1D* histNSHG_Sub[gMaxBoard][gMaxChannels];            // high gain all triggers BG sub
0702     TH2D* histNSLGHG[gMaxBoard][gMaxChannels];              // LG vs HG correlation
0703     TH2D* histNSHGLG[gMaxBoard][gMaxChannels];              // HG vs LG correlation
0704     TSpectrum* spectrumFitter[gMaxBoard][gMaxChannels];     // HG spectrum fittergHeader
0705     TF1* fitMultGaussNSHG[gMaxBoard][gMaxChannels];         // fit HG multi gauss
0706     TF1* fitMultGaussNSHGPreset[gMaxBoard][gMaxChannels];   // fit HG multi gauss preset values for peaks
0707     //***************************************************
0708     // mapped channels - layer & read-out board channels
0709     //***************************************************
0710     // channel counting in layer starts with 1, easier acces make array go to 8
0711     TH1D* histNSLG_mapped[gMaxLayers][9];                   // low gain all triggers
0712     TH1D* histNSHG_mapped[gMaxLayers][9];                   // high gain all triggers
0713     TH1D* histNSHG_BG_mapped[gMaxLayers][9];                // high gain all triggers background mapped
0714     TH1D* histNSHG_Sub_mapped[gMaxLayers][9];               // high gain all triggers background sub
0715     TF1* fitMultGaussNSHG_mapped[gMaxLayers][9];            // fit HG multi gauss
0716     TF1* fitMultGaussNSHGPreset_mapped[gMaxLayers][9];            // fit HG multi gauss
0717     //***************************************************
0718     // mapped channels - layer & read-out board channels - rebinned
0719     //***************************************************
0720     // channel counting in layer starts with 1, easier acces make array go to 8
0721     TH1D* histNSHG_mappedReb[gMaxLayers][9];                   // high gain all triggers
0722      
0723     //***************************************************
0724     // map of triggers above threshold  
0725     //***************************************************
0726     TH3D* histNS3DMap = new TH3D("h_map_NS_z_x_y", "; layer; col; row", 14, -0.5, 13.5, 4, -0.5, 3.5, 2, -0.5, 1.5);
0727     TH2D* histNS2DMap = new TH2D("h_map_NS_z_channel", "; channel; layer", 8, 0.5, 8.5, 14, -0.5, 13.5);
0728     TH1D* histNS1DMap = new TH1D("h_map_NS_channel", "; channel", 8, 0.5, 8.5);
0729     //***************************************************
0730     // Init histograms
0731     //***************************************************    
0732     for (Int_t j = 0; j < gMaxBoard; j++){
0733       for (Int_t i = 0; i < gMaxChannels; i++){
0734         histNSHG[j][i]    = new TH1D(Form("h_NS_HG_B%d_C%02d",j,i),";HG (adc); counts",4201,-200,4001);
0735         histNSLG[j][i]    = new TH1D(Form("h_NS_LG_B%d_C%02d",j,i),";LG (adc); counts",4201,-200,4001);
0736         histNSLGHG[j][i]  = new TH2D(Form("h_NS_LGHG_B%d_C%02d",j,i),";LG (adc); HG (adc)",420,-200,4001,420,-200,4001);
0737         histNSHGLG[j][i]  = new TH2D(Form("h_NS_HGLG_B%d_C%02d",j,i),";HG (adc); LG (adc)",420,-200,4001,420,-200,4001);
0738         spectrumFitter[j][i]  = new TSpectrum();
0739         histNSHG_BG[j][i] = nullptr;
0740       }
0741     }
0742     for (Int_t l = 0; l < gMaxLayers; l++){
0743       for (Int_t c = 0; c < 9; c++){
0744         histNSLG_mapped[l][c]             = nullptr;
0745         histNSHG_mapped[l][c]             = nullptr; 
0746         histNSHG_BG_mapped[l][c]          = nullptr; 
0747         histNSHG_mappedReb[l][c]          = nullptr;
0748         fitMultGaussNSHG_mapped[l][c]     = nullptr;
0749         fitMultGaussNSHGPreset_mapped[l][c]     = nullptr;
0750       }
0751     }
0752 
0753     // rebin array
0754     Double_t binningADC[3000];
0755     for (Int_t i = 0; i < 1200; i++) binningADC[i]                      = -200+i;
0756     for (Int_t i = 0; i < 500; i++) binningADC[i+1200]                  = 1000+i*2;   
0757     for (Int_t i = 0; i < 400; i++) binningADC[i+1200+500]              = 2000+i*5;   
0758     for (Int_t i = 0; i < 200; i++) binningADC[i+1200+500+400]          = 4000+i*10;   
0759     for (Int_t i = 0; i < 280; i++) binningADC[i+1200+500+400+200]      = 6000+i*50;   
0760     for (Int_t i = 0; i < 201; i++) binningADC[i+1200+500+400+200+280]  = 20000+i*100;   
0761     
0762     if (verbosity > 2)for (Int_t i = 0; i < 2781; i++) std::cout << binningADC[i] << "," ;
0763     if (verbosity > 2)std::cout<< std::endl; 
0764     
0765     // 1D channel representation of fit values, x axis scales as 10x layer count + channel within one assembley, 
0766     // - layer 3 channel 3: 33
0767     // - layer 0 channel 2: 2
0768     TH1D* hist1DnSPEPeaks_HG          = new TH1D("hist1DnSPEPeaks_HG_channels", "; 10x layer + board channel; # SPE peaks", 10*tempLayerBinning+1, -0.5, 10*tempLayerBinning+0.5 );
0769     TH1D* hist1DAvDiffSPEPeaks_HG      = new TH1D("hist1DAvDiffSPEPeaks_HG_channels", "; 10x layer + board channel; #mu(#Delta_{SPE}) (HG ADC)", 10*tempLayerBinning+1, -0.5, 10*tempLayerBinning+0.5 );
0770     TH1D* hist1DAvDiffSPEPeaksFit_HG   = new TH1D("hist1DAvDiffSPEPeaksFit_HG_channels", "; 10x layer + board channel; #mu(#Delta_{SPE,fit}) (HG ADC)", 10*tempLayerBinning+1, -0.5, 10*tempLayerBinning+0.5 );
0771 
0772     // 1D channel representation of fit values, x axis scales as CAEN board channelns 64x CAEN board # + CAEN channel
0773     TH1D* hist1DCAEN_nSPEPeaks_HG           = new TH1D("hist1DCAEN_nSPEPeaks_HG_channels", "; 64x CAEN board + CAEN channel; # SPE peaks", 64*gMaxBoard+1, -0.5, 64*gMaxBoard+0.5 );
0774     TH1D* hist1DCAEN_AvDiffSPEPeaks_HG      = new TH1D("hist1DCAEN_AvDiffSPEPeaks_HG_channels", "; 64x CAEN board + CAEN channel; #mu(#Delta_{SPE}) (HG ADC)", 64*gMaxBoard+1, -0.5, 64*gMaxBoard+0.5 );
0775     TH1D* hist1DCAEN_AvDiffSPEPeaksFit_HG   = new TH1D("hist1DCAEN_AvDiffSPEPeaksFit_HG_channels", "; 64x CAEN board + CAEN channel; #mu(#Delta_{SPE,fit}) (HG ADC)", 64*gMaxBoard+1, -0.5, 64*gMaxBoard+0.5 );
0776         
0777     // 2D representation of fit values
0778     TH2D* hist2DnSPEPeaks_HG          = new TH2D("hist2DnSPEPeaks_HG_z_channel", "; channel; layer; # SPE (HG ADC)", 8, 0.5, 8.5, 14, -0.5, 13.5);
0779     TH2D* hist2DAvDiffSPEPeaks_HG     = new TH2D("hist2DAvDiffSPEPeaks_HG_z_channel", "; channel; layer; #mu(#Delta_{SPE}) (HG ADC)", 8, 0.5, 8.5, 14, -0.5, 13.5);
0780     TH2D* hist2DAvDiffSPEPeaksFits_HG = new TH2D("hist2DAvDiffSPEPeaksFits_HG_z_channel", "; channel; layer; #mu(#Delta_{SPE}) (HG ADC)", 8, 0.5, 8.5, 14, -0.5, 13.5);
0781     
0782     Int_t nSPE[gMaxLayers][9]                = {{0}};
0783     Double_t avDiffSPEPeaks[gMaxLayers][9]      = {{0.}};
0784     Double_t avDiffSPEPeaksFits[gMaxLayers][9]  = {{0.}};
0785   
0786     // ********************************************************************************************************
0787     // Second loop over full tree to obtain noise subtracted histograms
0788     // ********************************************************************************************************
0789     TRandom3* rand    = new TRandom3();
0790     nEventsProcessed  = 0;
0791     for (Long64_t i=startEvent; i<nEntriesTree;i++) {
0792         // load current event
0793         tt_event->GetEntry(i);
0794         nEventsProcessed++;
0795 
0796         Int_t nChNoNoise = 0;
0797         // processing progress info
0798         if(i>0 && nEntriesTree>100 && i%(nEntriesTree/(20))==0) std::cout << "//processed " << 100*(i)/nEntriesTree << "%"  << std::endl;
0799         if(verbosity>1){
0800           std::cout << "***********************************************************************************************************" << std::endl;
0801           std::cout << "event " << i << std::endl;
0802           std::cout << "***********************************************************************************************************" << std::endl;
0803         }
0804         // temporary mapping of channels after noise subtractions
0805         Int_t signal[2][64][4] = {{{0}}};         // [board][channel][HG, LG, trigg HG,comb HG&LG]
0806         Float_t signalNS[2][64][4] = {{{0}}};       // [board][channel][HG, LG, trigg HG,comb HG&LG]
0807         Double_t summed[2] = {0};                 // total channel sum
0808         // readEvent
0809         if (verbosity > 1)std::cout << "------------- Event -------------------------" << std::endl;
0810         if (verbosity > 1)std::cout << gTrID << "\t" << gTRtimeStamp << std::endl;
0811         if (i == startEvent) tstapMin = gTRtimeStamp;
0812         tstapMax = gTRtimeStamp;
0813         Double_t tCurr= (tstapMax-tstapMin)/1e6/60; // elapsed time in min
0814         
0815         for(Int_t ch=0; ch<gMaxChannels && ch < nChmapped; ch++){
0816           Double_t meanNHG = 0;
0817           if (mean[0][gBoard[ch]][gChannel[ch]] != -1)
0818             meanNHG = mean[0][gBoard[ch]][gChannel[ch]];
0819           Double_t meanNLG = 0;
0820           if (mean[2][gBoard[ch]][gChannel[ch]] != -1)
0821             meanNLG = mean[2][gBoard[ch]][gChannel[ch]];
0822 
0823           if (i == 0){
0824             std::cout << ch << "\t means: " <<  meanNHG << "\t" << meanNLG << std::endl;
0825           }
0826           
0827           signal[gBoard[ch]][gChannel[ch]][0] = gHG[ch];
0828           signal[gBoard[ch]][gChannel[ch]][1] = gLG[ch];
0829           
0830           signalNS[gBoard[ch]][gChannel[ch]][0] = gHG[ch]-meanNHG;
0831           signalNS[gBoard[ch]][gChannel[ch]][1] = gLG[ch]-meanNLG;
0832           // set combined HG & LG signal
0833           if (signalNS[gBoard[ch]][gChannel[ch]][0] < 3800)
0834             signalNS[gBoard[ch]][gChannel[ch]][3] = signalNS[gBoard[ch]][gChannel[ch]][0];
0835           else 
0836             signalNS[gBoard[ch]][gChannel[ch]][3] = signalNS[gBoard[ch]][gChannel[ch]][1] *cslope[0][gBoard[ch]][gChannel[ch]] + rand->Rndm()*cslope[0][gBoard[ch]][gChannel[ch]];
0837           Double_t chThres = scaledThr-meanNHG;
0838           if (gHG[ch] > chThres)  // fill potential trigger signals
0839             signalNS[gBoard[ch]][gChannel[ch]][2] = gHG[ch]-chThres;
0840           // fill histos
0841           histNSHG[gBoard[ch]][gChannel[ch]]->Fill(signalNS[gBoard[ch]][gChannel[ch]][0]); 
0842           histNSLG[gBoard[ch]][gChannel[ch]]->Fill(signalNS[gBoard[ch]][gChannel[ch]][1]); 
0843           histNSLGHG[gBoard[ch]][gChannel[ch]]->Fill(signalNS[gBoard[ch]][gChannel[ch]][1], signalNS[gBoard[ch]][gChannel[ch]][0]); 
0844           histNSHGLG[gBoard[ch]][gChannel[ch]]->Fill(signalNS[gBoard[ch]][gChannel[ch]][0], signalNS[gBoard[ch]][gChannel[ch]][1]); 
0845           if (signalNS[gBoard[ch]][gChannel[ch]][0] > (adcThreshold-meanNHG)){
0846             nChNoNoise++;
0847             if (verbosity > 1) std::cout << "not noise" << std::endl;
0848             histNS3DMap->Fill(mapping[gChannel[ch]][0],mapping[gChannel[ch]][3],mapping[gChannel[ch]][2]);
0849             histNS2DMap->Fill(mapping[gChannel[ch]][1],mapping[gChannel[ch]][0]);
0850             histNS1DMap->Fill(mapping[gChannel[ch]][1]);
0851           }
0852           // calculate summed signal
0853           summed[gBoard[ch]] = summed[gBoard[ch]]+signalNS[gBoard[ch]][gChannel[ch]][3]; 
0854         }
0855         // fill summed signal
0856         if (verbosity > 1)std::cout << "Channels above noise: "<< nChNoNoise << std::endl;        
0857     }
0858     std::cout << "Total events processed:" << nEventsProcessed << std::endl;
0859 
0860     // ********************************************************************************************************
0861     // Plotting of single channels after noise subtraction
0862     // -> remapping of CAEN board+channels to physicsal layers and readout-channels
0863     // ********************************************************************************************************    
0864     for (Int_t j = 0; j < gMaxBoard; j++){
0865       for (Int_t i = 0; i < gMaxChannels; i++){
0866         Int_t chMap     = j*64 + i;
0867         Int_t chBoard   = mapping[chMap][1];
0868         Int_t layer     = mapping[chMap][0];
0869         if (! ((histNSHG[j][i] && histNSHG[j][i]->GetEntries() > 0)|| (histNSLG[j][i] && histNSLG[j][i]->GetEntries() > 0)) ){
0870           fitMultGaussNSHG[j][i] = nullptr;
0871           fitMultGaussNSHG_mapped[layer][chBoard] = nullptr;
0872           continue;
0873         }
0874         if (verbosity > 0)std::cout << j << "\t" << i << "\t" << chMap << "\t L: " << layer  << "\t C:" <<  chBoard << std::endl;
0875         
0876         if (currentRunInfo.vop < 3.5){
0877           histNSHG_BG[j][i]                   = (TH1D*)spectrumFitter[j][i]->Background(histNSHG[j][i], 14, "smoothing3");
0878         } else {
0879            histNSHG_BG[j][i]                   = (TH1D*)spectrumFitter[j][i]->Background(histNSHG[j][i], 18, "smoothing3");
0880         }
0881         histNSHG_BG_mapped[layer][chBoard]  = (TH1D*)histNSHG_BG[j][i]->Clone(Form("h_NS_HG_BG_mapped_L%d_C%02d",layer,chBoard));
0882         histNSHG_Sub[j][i]                  = (TH1D*)histNSHG[j][i]->Clone(Form("h_NS_HG_Sub_B%d_C%02d",j,i));
0883         histNSHG_Sub[j][i]->Sumw2();
0884         histNSHG_Sub[j][i]->Add(histNSHG_BG[j][i], -1);
0885         histNSHG_Sub_mapped[layer][chBoard]  = (TH1D*)histNSHG_Sub[j][i]->Clone(Form("h_NS_HG_Sub_mapped_L%d_C%02d",layer,chBoard));
0886         
0887         if (currentRunInfo.vop < 3.5){
0888           nSPE[j][i] = spectrumFitter[j][i]->Search(histNSHG_Sub_mapped[layer][chBoard], 5., "", 0.05); // Search for peaks
0889         } else {
0890           nSPE[j][i] = spectrumFitter[j][i]->Search(histNSHG_Sub_mapped[layer][chBoard], 7., "", 0.05); // Search for peaks
0891         }
0892         if (verbosity > 0) std::cout << "found " << nSPE[j][i] << "\t different peaks in spectrum" << std::endl;
0893         
0894         if (nSPE[j][i] > 2){
0895           nPeaksMultGauss = 0;
0896           Double_t *xpeaks;
0897           Double_t par[3000];
0898           Double_t par2[3000];
0899           xpeaks = spectrumFitter[j][i]->GetPositionX();
0900           
0901           // invert order of peaks
0902           std::sort(xpeaks, xpeaks + nSPE[j][i], [&](Double_t a, Double_t b) {
0903               return a > b;
0904           });
0905           
0906           std::cout << "Differences between consecutive entries:" << std::endl;
0907           Double_t sum_diff = 0.0;
0908           for (size_t f = 1; f < (size_t)nSPE[j][i]-1; ++f) {
0909             Double_t diff = xpeaks[f - 1] - xpeaks[f];
0910             std::cout << diff << " ";
0911             sum_diff += diff;
0912           }
0913           std::cout << std::endl;
0914           
0915           avDiffSPEPeaks[j][i] = sum_diff / (nSPE[j][i] - 2);
0916           std::cout << "Average of differences: " << avDiffSPEPeaks[j][i] << std::endl;
0917           par[0]  = 0;
0918           par[1]  = 0;
0919           
0920           for (Int_t p = 0; p < (Int_t)nSPE[j][i]; p++){
0921             Double_t xp   = xpeaks[p];
0922 
0923             Int_t bin     = histNSHG_Sub[j][i]->GetXaxis()->FindBin(xp);
0924             Double_t yp   = histNSHG_Sub[j][i]->GetBinContent(bin);
0925 
0926             if (verbosity > 1) std::cout << p << "\t" << xp << "\t" << bin << "\t"<< yp << std::endl;            
0927             par[3 * nPeaksMultGauss + 2] = yp;     // "height"
0928             par[3 * nPeaksMultGauss + 3] = xp; // "mean"
0929             // par[3 * nPeaksMultGauss + 4] = 5;  // "sigma"
0930             par[3 * nPeaksMultGauss + 4] = sigma[0][j][i];  // "sigma"
0931             if (verbosity > 1) std::cout<<xp<<std::endl;
0932             nPeaksMultGauss++;
0933           }
0934           
0935           fitMultGaussNSHG[j][i] = new TF1(Form("fitMultGauss_NS_HG_Sub_B%d_C%02d",j,i), multGauss, 0, 4096, 3 * nPeaksMultGauss+2);
0936           fitMultGaussNSHG[j][i]->SetParameters(par);
0937           fitMultGaussNSHG[j][i]->SetNpx(1000);
0938           fitMultGaussNSHGPreset[j][i] = new TF1(Form("fitMultGauss_NS_HG_SubPre_B%d_C%02d",j,i), multGauss, 0, 4096, 3 * nPeaksMultGauss+2);
0939           fitMultGaussNSHGPreset[j][i]->SetParameters(par);
0940           for (Int_t p = 0; p < (Int_t)nSPE[j][i]; p++){
0941               fitMultGaussNSHG[j][i]->SetParLimits( 3*p+3, xpeaks[p] - 20, xpeaks[p] + 20);
0942               fitMultGaussNSHGPreset[j][i]->FixParameter(3*p+3, fitMultGaussNSHGPreset[j][i]->GetParameter(3*p+3));
0943               
0944           }
0945           
0946           fitMultGaussNSHGPreset[j][i]->SetNpx(1000);
0947           if (verbosity > 1) {
0948             for (Int_t n = 0; n < fitMultGaussNSHG[j][i]->GetNpar(); n++){
0949                 std::cout << n << "\t" << fitMultGaussNSHG[j][i]->GetParameter(n) << std::endl;
0950             }
0951           }
0952           histNSHG_Sub[j][i]->Fit(fitMultGaussNSHG[j][i],"QRMNE0");
0953           histNSHG_Sub[j][i]->Fit(fitMultGaussNSHGPreset[j][i],"QRMNE0");
0954           Double_t sum_diff_Fit     = 0.0;
0955           Double_t xPeaksFits[3000] = {0};
0956           for (Int_t p = 0; p < (Int_t)nSPE[j][i]-1; p++){
0957             std::cout << p << "\t mu=" <<  fitMultGaussNSHG[j][i]->GetParameter(3*p + 3) << "\t sigma=" << fitMultGaussNSHG[j][i]->GetParameter(3*p + 4);
0958             xPeaksFits[p] = fitMultGaussNSHG[j][i]->GetParameter(3*p + 3);
0959             if (p > 0) {
0960                 Double_t diffTemp = xPeaksFits[p-1] - xPeaksFits[p];
0961                 std::cout << "\t diff = " << diffTemp;
0962                 sum_diff_Fit += diffTemp;
0963             }
0964             std::cout << std::endl;
0965           }
0966           avDiffSPEPeaksFits[j][i] = sum_diff_Fit / (nSPE[j][i] - 2);
0967           std::cout << "Average of differences fit: " << avDiffSPEPeaksFits[j][i] << std::endl;
0968         } else {
0969           std::cout << "no peaks from different pixels found" << std::endl;
0970           fitMultGaussNSHG[j][i] = nullptr;
0971           fitMultGaussNSHGPreset[j][i] = nullptr;
0972         }
0973 
0974         // ****************************************************************************
0975         // fill Monitoring histogramms
0976         // ****************************************************************************
0977         hist2DnSPEPeaks_HG->Fill(chBoard, layer,nSPE[j][i]);
0978         hist2DAvDiffSPEPeaks_HG->Fill(chBoard, layer,avDiffSPEPeaks[j][i]);
0979         hist2DAvDiffSPEPeaksFits_HG->Fill(chBoard, layer,avDiffSPEPeaks[j][i]);
0980         
0981         /// get bin to fill for 1D representation of channels for fit values
0982         Int_t channelBin1D = hist1DnSPEPeaks_HG->FindBin(layer*10+chBoard);          
0983         hist1DnSPEPeaks_HG->SetBinContent(channelBin1D, nSPE[j][i]);
0984         hist1DAvDiffSPEPeaks_HG->SetBinContent(channelBin1D, avDiffSPEPeaks[j][i]);
0985         hist1DAvDiffSPEPeaksFit_HG->SetBinContent(channelBin1D, avDiffSPEPeaksFits[j][i]);
0986         
0987         hist1DCAEN_nSPEPeaks_HG->SetBinContent(chMap+1, nSPE[j][i]);
0988         hist1DCAEN_AvDiffSPEPeaks_HG->SetBinContent(chMap+1, avDiffSPEPeaks[j][i]);
0989         hist1DCAEN_AvDiffSPEPeaksFit_HG->SetBinContent(chMap+1, avDiffSPEPeaksFits[j][i]);
0990         
0991         // ****************************************************************************
0992         // create layer mapped hists/fits
0993         // ****************************************************************************
0994         histNSHG_mapped[layer][chBoard] = (TH1D*)histNSHG[j][i]->Clone(Form("h_NS_HG_mapped_L%d_C%02d",layer,chBoard));
0995         histNSLG_mapped[layer][chBoard] = (TH1D*)histNSLG[j][i]->Clone(Form("h_NS_LG_mapped_L%d_C%02d",layer,chBoard));
0996         if (fitMultGaussNSHG[j][i]){
0997           fitMultGaussNSHG_mapped[layer][chBoard] = (TF1*)fitMultGaussNSHG[j][i]->Clone(Form("fitMultiGauss_NS_HG_mapped_L%d_C%02d",layer,chBoard));
0998           fitMultGaussNSHGPreset_mapped[layer][chBoard] = (TF1*)fitMultGaussNSHG[j][i]->Clone(Form("fitMultiGauss_NS_HGPre_mapped_L%d_C%02d",layer,chBoard));
0999         }
1000         histNSHG_mapped[layer][chBoard]->Sumw2();
1001         histNSHG_mappedReb[layer][chBoard] = (TH1D*)histNSHG_mapped[layer][chBoard]->Rebin(2100,Form("%sReb",(histNSHG_mapped[layer][chBoard]->GetName())), binningADC);
1002         histNSHG_mappedReb[layer][chBoard]->Scale(1,"width");
1003       
1004         if (bDetPlot){
1005           // ********************************************************
1006           // Plot different triggers together
1007           // ********************{************************************
1008           PlotOverlayDiffTriggers( canvas1DDiffTrigg, histNSHG[j][i], nullptr, nullptr, NULL,
1009                                 -100, 1100, Form("%s/HG_NS_DiffTriggers", outputDirPlotsDet.Data()),
1010                                 j, i, layer, chBoard, currentRunInfo, 0.04);
1011           PlotOverlayDiffTriggers( canvas1DDiffTrigg, histNSLG[j][i], nullptr, nullptr, NULL,
1012                                 -100, 500, Form("%s/LG_NS_DiffTriggers", outputDirPlotsDet.Data()),
1013                                 j, i, layer, chBoard, currentRunInfo, 0.04);
1014           PlotOverlaySinglePhoton( canvas1DDiffTrigg, histNSHG[j][i], histNSHG_BG[j][i], histNSHG_Sub[j][i],  fitMultGaussNSHG[j][i],
1015                                 -100, 1100, Form("%s/HG_NS_WithSpectrumBG_DiffTriggers", outputDirPlotsDet.Data()),
1016                                 j, i, layer, chBoard, currentRunInfo, 0.04);
1017           PlotOverlaySinglePhoton( canvas1DDiffTrigg, histNSHG[j][i], histNSHG_BG[j][i], histNSHG_Sub[j][i],  fitMultGaussNSHGPreset[j][i],
1018                                 -100, 1100, Form("%s/HG_NS_WithSpectrumBGPreset_DiffTriggers", outputDirPlotsDet.Data()),
1019                                 j, i, layer, chBoard, currentRunInfo, 0.04);
1020         }
1021       }
1022     }
1023 
1024     PlotSimpleMultiLayer2D( canvas2DCorr, hist2DNoiseMean_HG, maxActiveLayer, maxActiveRBCh, textSizeRel, Form("%s/HG_Noise_Mean.pdf", outputDirPlots.Data()), currentRunInfo);
1025     PlotSimpleMultiLayer2D( canvas2DCorr, hist2DNoiseSigma_HG, maxActiveLayer, maxActiveRBCh, textSizeRel,Form("%s/HG_Noise_Sigma.pdf", outputDirPlots.Data()), currentRunInfo);
1026     PlotSimpleMultiLayer2D( canvas2DCorr, hist2DNoiseMean_LG, maxActiveLayer, maxActiveRBCh, textSizeRel,Form("%s/LG_Noise_Mean.pdf", outputDirPlots.Data()), currentRunInfo);
1027     PlotSimpleMultiLayer2D( canvas2DCorr, hist2DNoiseSigma_LG, maxActiveLayer, maxActiveRBCh, textSizeRel, Form("%s/LG_Noise_Sigma.pdf", outputDirPlots.Data()), currentRunInfo);
1028     PlotSimpleMultiLayer2D( canvas2DCorr, hist2DLGHG_slope, maxActiveLayer, maxActiveRBCh, textSizeRel, Form("%s/LGHG_Slope.pdf", outputDirPlots.Data()), currentRunInfo);
1029     PlotSimpleMultiLayer2D( canvas2DCorr, hist2DnSPEPeaks_HG, maxActiveLayer, maxActiveRBCh, textSizeRel, Form("%s/NPeaksVisible.pdf", outputDirPlots.Data()), currentRunInfo, kTRUE);
1030     PlotSimpleMultiLayer2D( canvas2DCorr, hist2DAvDiffSPEPeaks_HG, maxActiveLayer, maxActiveRBCh, textSizeRel, Form("%s/AvDiffSPEPeaks.pdf", outputDirPlots.Data()), currentRunInfo, kTRUE);
1031     PlotSimpleMultiLayer2D( canvas2DCorr, hist2DAvDiffSPEPeaksFits_HG, maxActiveLayer, maxActiveRBCh, textSizeRel, Form("%s/AvDiffSPEPeaksFits.pdf", outputDirPlots.Data()), currentRunInfo, kTRUE);
1032     
1033     
1034     for (Int_t l = 0; l < gMaxLayers; l++){
1035       if (!lActive[l]) continue; 
1036       PlotDiffTriggersFullLayer (canvas8Panel,pad8Panel, topRCornerX, topRCornerY, relSize8P, textSizePixel, 
1037                                 histNSHG_mapped[l], nullptr, nullptr, nullptr, 
1038                                 -80, 1000, 1.2, l , Form("%s/TriggerOverlay_HG_NS_Zoomed_Layer%02d.pdf" ,outputDirPlots.Data(), l), currentRunInfo);
1039       PlotSPEFullLayer (canvas8Panel,pad8Panel, topRCornerX, topRCornerY, relSize8P, textSizePixel, 
1040                                 histNSHG_mapped[l], histNSHG_Sub_mapped[l], histNSHG_BG_mapped[l], fitMultGaussNSHG_mapped[l],
1041                                 -80, 500, 1.2, l , Form("%s/SPEOverlay_HG_NS_Zoomed_Layer%02d.pdf" ,outputDirPlots.Data(), l), currentRunInfo);
1042       PlotSPEFullLayer (canvas8Panel,pad8Panel, topRCornerX, topRCornerY, relSize8P, textSizePixel, 
1043                                 histNSHG_mapped[l], histNSHG_Sub_mapped[l], histNSHG_BG_mapped[l], fitMultGaussNSHGPreset_mapped[l],
1044                                 -80, 500, 1.2, l , Form("%s/SPEOverlay_HG_NS_Preset_Zoomed_Layer%02d.pdf" ,outputDirPlots.Data(), l), currentRunInfo);
1045       PlotDiffTriggersFullLayer (canvas8Panel,pad8Panel, topRCornerX, topRCornerY, relSize8P, textSizePixel, 
1046                                 histNSLG_mapped[l],  nullptr, nullptr, nullptr, 
1047                                 -80, 500, 1.2, l , Form("%s/TriggerOverlay_LG_NS_Zoomed_Layer%02d.pdf" ,outputDirPlots.Data(), l), currentRunInfo);
1048     }
1049     
1050     // ********************************************************************************************************
1051     // Overlay noise in same layer
1052     // ********************************************************************************************************
1053     for (Int_t l = 0; l < gMaxLayers; l++){
1054       if (!lActive[l]) continue;      
1055 //       for (Int_t c = 1; c < 9; c++){
1056 //         
1057 //         PlotMIPSingle (canvas1DNoise, histNSHG_mappedReb[l][c],nullptr, chisqr[l][c], ndf[l][c], maxLandG[l][c], fwhmLandG[l][c],
1058 //                         l, c, Form("%s/HG_TriggMipWithFit", outputDirPlotsDet.Data()), currentRunInfo, 0.04);
1059 //       }
1060       PlotChannelOverlaySameLayer( canvas1DDiffTrigg, histNSHG_mappedReb[l], 1, 9,
1061                                   -100, 1000, 1./5, Form("%s/HGNS", outputDirPlots.Data()), l, currentRunInfo, 0.04,"hist");
1062       PlotChannelOverlaySameLayer( canvas1DDiffTrigg, histNSLG_mapped[l], 1, 9,
1063                                   -100, 500, 1, Form("%s/LGNS", outputDirPlots.Data()), l, currentRunInfo, 0.04);        
1064     }
1065     
1066     // ********************************************************************************************************
1067     // Overlay noise for same readout channel
1068     // ********************************************************************************************************
1069     for (Int_t c = 1; c < 9; c++){      
1070       PlotChannelOverlaySameReadoutChannel( canvas1DDiffTrigg, histNSHG_mappedReb, lActive, 0, gMaxLayers, 
1071                                             -100,1000, 1./5, Form("%s/HGNS", outputDirPlots.Data()), c, currentRunInfo,0.04,"hist");    
1072     }
1073 
1074     // *****************************************************************
1075     // Create graphs per board and channels with calib values MIP values
1076     // *****************************************************************
1077     // only create data point if more than 2 peaks were visible in SPE in graphs
1078     TGraphErrors* graphnSPEPeaks_HG           = CreateGraphFromHistAndCleanup(hist1DnSPEPeaks_HG, "graphnSPEPeaks_HG_channels", 3, kTRUE);
1079     TGraphErrors* graphAvDiffSPEPeaks_HG      = CreateGraphFromHistAndCleanup(hist1DAvDiffSPEPeaks_HG, "graphAvDiffSPEPeaks_HG_channels");
1080     TGraphErrors* graphAvDiffSPEPeaksFit_HG   = CreateGraphFromHistAndCleanup(hist1DAvDiffSPEPeaksFit_HG, "graphAvDiffSPEPeaksFit_HG_channels");
1081 
1082     TGraphErrors* graphCAEN_nSPEPeaks_HG           = CreateGraphFromHistAndCleanup(hist1DCAEN_nSPEPeaks_HG, "graphCAEN_nSPEPeaks_HG_channels", 3, kTRUE);
1083     TGraphErrors* graphCAEN_AvDiffSPEPeaks_HG      = CreateGraphFromHistAndCleanup(hist1DCAEN_AvDiffSPEPeaks_HG, "graphCAEN_AvDiffSPEPeaks_HG_channels");
1084     TGraphErrors* graphCAEN_AvDiffSPEPeaksFit_HG   = CreateGraphFromHistAndCleanup(hist1DCAEN_AvDiffSPEPeaksFit_HG, "graphCAEN_AvDiffSPEPeaksFit_HG_channels");
1085 
1086     TObjString* stringRunInfo = new TObjString;
1087     stringRunInfo->SetString(GetStringFromRunInfo(currentRunInfo,4));
1088     
1089     // ********************************************************************************************************
1090     // write output to single file
1091     // ********************************************************************************************************
1092     TFile* fileOutput = new TFile(Form("%s/output_basics.root",outputDir.Data()),"RECREATE");
1093     // ******************************************************
1094     // create folders
1095     // ******************************************************
1096     fileOutput->mkdir("RawData");
1097     fileOutput->mkdir("RawDataMapped");
1098     fileOutput->mkdir("NoiseSubtracted");
1099     fileOutput->mkdir("NoiseSubtractedMapped");
1100     fileOutput->mkdir("CorrelationLGHG");
1101     fileOutput->mkdir("CorrelationHGLG");
1102     // ******************************************************
1103     // run over all CAEN boards and channels
1104     // ******************************************************
1105     for (Int_t j = 0; j < gMaxBoard; j++){
1106       for (Int_t i = 0; i < gMaxChannels; i++){
1107         // raw data
1108         fileOutput->cd("RawData");
1109         WriteOnlyIfFilled(histHG[j][i]);
1110         WriteOnlyIfFilled(hist_T_HG[j][i]);
1111         WriteOnlyIfFilled(histLG[j][i]);
1112         // noise subtracted
1113         fileOutput->cd("NoiseSubtracted");
1114         WriteOnlyIfFilled(histNSHG[j][i]);
1115         WriteOnlyIfFilled(histNSLG[j][i]);
1116         if (fitMultGaussNSHG[j][i]) fitMultGaussNSHG[j][i]->Write();
1117         // correlation histograms
1118         fileOutput->cd("CorrelationLGHG");
1119         WriteOnlyIfFilled(histLGHG[j][i]);
1120         WriteOnlyIfFilled(histNSLGHG[j][i]);
1121         if (fitLGHGCorr[j][i]) fitLGHGCorr[j][i]->Write();
1122         fileOutput->cd("CorrelationHGLG");
1123         WriteOnlyIfFilled(histHGLG[j][i]);
1124         WriteOnlyIfFilled(histNSHGLG[j][i]);
1125         if (fitHGLGCorr[j][i]) fitHGLGCorr[j][i]->Write();
1126       }
1127       fileOutput->cd();
1128     }
1129     std::cout << __LINE__ << std::endl;
1130     // *******************************************************
1131     // run over all physical layers and readout board channels
1132     // *******************************************************
1133     for (Int_t l = 0; l  < gMaxLayers; l++){
1134       for (Int_t c = 1; c < 9; c++){
1135         fileOutput->cd("RawDataMapped");
1136         WriteOnlyIfFilled(histHG_mapped[l][c]);
1137         if (fitGausHG_BG_mapped[l][c]) fitGausHG_BG_mapped[l][c]->Write();
1138         WriteOnlyIfFilled(histLG_mapped[l][c]);
1139         if (fitGausLG_BG_mapped[l][c]) fitGausLG_BG_mapped[l][c]->Write();
1140         fileOutput->cd("NoiseSubtractedMapped");
1141         WriteOnlyIfFilled(histNSHG_mapped[l][c]);
1142         WriteOnlyIfFilled(histNSHG_mappedReb[l][c]);
1143         WriteOnlyIfFilled(histNSLG_mapped[l][c]);
1144       }
1145     }
1146     // *******************************************************
1147     // write summary histograms
1148     // *******************************************************
1149     fileOutput->cd();
1150     stringRunInfo->Write();
1151     histNChAboveNoise->Write();
1152     hist2DNoiseMean_HG->Write();
1153     histNoiseMean_HG->Write();
1154     hist2DNoiseSigma_HG->Write();
1155     histNoiseSigma_HG->Write();
1156     hist2DNoiseMean_LG->Write();
1157     histNoiseMean_LG->Write();
1158     hist2DNoiseSigma_LG->Write();
1159     histNoiseSigma_LG->Write();
1160     hist2DLGHG_slope->Write();
1161     histLGHG_slope->Write();
1162     hist2DLGHG_offset->Write();
1163     
1164     hist1DNoiseSigma_HG->Write();
1165     hist1DNoiseMean_HG->Write();
1166     hist1DNoiseSigma_LG->Write();
1167     hist1DNoiseMean_LG->Write();
1168     hist1DLGHG_slope->Write();
1169     hist1DLGHG_offset->Write();
1170     hist1DHGLG_slope->Write();
1171     hist1DHGLG_offset->Write();
1172     hist1DCAEN_NoiseSigma_HG->Write();
1173     hist1DCAEN_NoiseMean_HG->Write();
1174     hist1DCAEN_NoiseSigma_LG->Write();
1175     hist1DCAEN_NoiseMean_LG->Write();
1176     hist1DCAEN_LGHG_slope->Write();
1177     hist1DCAEN_LGHG_offset->Write();
1178     hist1DCAEN_HGLG_slope->Write();
1179     hist1DCAEN_HGLG_offset->Write();
1180 
1181     gNoiseMeanHG->Write();
1182     gNoiseSigmaHG->Write();
1183     gNoiseMeanLG->Write();
1184     gNoiseSigmaLG->Write();
1185     gCorrLGHGSlope->Write();
1186     gCorrLGHGOffset->Write();
1187     gCorrHGLGSlope->Write();
1188     gCorrHGLGOffset->Write();
1189     gCAEN_NoiseMeanHG->Write();
1190     gCAEN_NoiseSigmaHG->Write();
1191     gCAEN_NoiseMeanLG->Write();
1192     gCAEN_NoiseSigmaLG->Write();
1193     gCAEN_CorrLGHGSlope->Write();
1194     gCAEN_CorrLGHGOffset->Write();
1195     gCAEN_CorrHGLGSlope->Write();
1196     gCAEN_CorrHGLGOffset->Write();
1197 
1198     hist1DnSPEPeaks_HG->Write();
1199     hist1DAvDiffSPEPeaks_HG->Write();
1200     hist1DAvDiffSPEPeaksFit_HG->Write();
1201     hist1DCAEN_nSPEPeaks_HG->Write();
1202     hist1DCAEN_AvDiffSPEPeaks_HG->Write();
1203     hist1DCAEN_AvDiffSPEPeaksFit_HG->Write();
1204 
1205     graphnSPEPeaks_HG->Write();
1206     graphAvDiffSPEPeaks_HG->Write();
1207     graphAvDiffSPEPeaksFit_HG->Write();
1208     graphCAEN_nSPEPeaks_HG->Write();
1209     graphCAEN_AvDiffSPEPeaks_HG->Write();
1210     graphCAEN_AvDiffSPEPeaksFit_HG->Write();
1211 
1212     hist2DnSPEPeaks_HG->Write();
1213     hist2DAvDiffSPEPeaks_HG->Write();
1214     hist2DAvDiffSPEPeaksFits_HG->Write();
1215 
1216     hist3DMap->Write();
1217     hist2DMap->Write();
1218     hist1DMap->Write();
1219     histNS3DMap->Write();
1220     histNS2DMap->Write();
1221     histNS1DMap->Write();
1222      // write output file
1223     fileOutput->Write();
1224     fileOutput->Close();
1225     // ********************************************************************************************************
1226     // DONE
1227     // ********************************************************************************************************
1228 
1229 }                                  
1230                                   
1231
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!