Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2026-06-06 07:54:26 

0001 #include "DataAnalysis.h"
0002 #include <vector>
0003 #include "TROOT.h"
0004 #include <cmath>  
0005 #include <algorithm>
0006 //#include <unistd.h> // Add for use on Mac OS -> Same goes for Analyse.cc
0007 #include "TF1.h"
0008 #include "TFitResult.h"
0009 #include "TFitResultPtr.h"
0010 #include "TH1D.h"
0011 #include "TH2D.h"
0012 #include "TProfile.h"
0013 #include "TChain.h"
0014 #include "CommonHelperFunctions.h"
0015 #include "TileSpectra.h"
0016 #include "MultiCanvas.h"
0017 #include "PlotHelper.h"
0018 #include "TRandom3.h"
0019 #include "TMath.h"
0020 #include "Math/Minimizer.h"
0021 #include "Math/MinimizerOptions.h"
0022 
0023 // ****************************************************************************
0024 // Checking and opening input and output files
0025 // ****************************************************************************
0026 bool DataAnalysis::CheckAndOpenIO(void){
0027   int matchingbranch;
0028 
0029   //Need to check first input to get the setup...I do not think it is necessary
0030     std::cout<<"Input name set to: '"<<RootInputName.Data() <<std::endl;  
0031   if(!RootInputName.IsNull()){
0032     //File exist?
0033     RootInput=new TFile(RootInputName.Data(),"READ");
0034     if(RootInput->IsZombie()){
0035       std::cout<<"Error opening '"<<RootInputName<<"', does the file exist?"<<std::endl;
0036       return false;
0037     }
0038 
0039     //Retrieve info, start with setup
0040     TsetupIn = (TTree*) RootInput->Get("Setup");
0041     if(!TsetupIn){
0042       std::cout<<"Could not retrieve the Setup tree, leaving"<<std::endl;
0043       return false;
0044     }
0045     setup=Setup::GetInstance();
0046     std::cout<<"Setup add "<<setup<<std::endl;
0047     //matchingbranch=TsetupIn->SetBranchAddress("setup",&setup);
0048     matchingbranch=TsetupIn->SetBranchAddress("setup",&rswptr);
0049     if(matchingbranch<0){
0050       std::cout<<"Error retrieving Setup info from the tree"<<std::endl;
0051       return false;
0052     }
0053     std::cout<<"Entries "<<TsetupIn->GetEntries()<<std::endl;
0054     TsetupIn->GetEntry(0);
0055     setup->Initialize(*rswptr);
0056     std::cout<<"Reading "<<RootInput->GetName()<<std::endl;
0057     std::cout<<"Setup Info "<<setup->IsInit()<<"  and  "<<setup->GetCellID(0,0)<<std::endl;
0058     //std::cout<<"Setup add now "<<setup<<std::endl;
0059 
0060     //Retrieve info, existing data?
0061     TdataIn = (TTree*) RootInput->Get("Data");
0062     if(!TdataIn){
0063       std::cout<<"Could not retrieve the Data tree, leaving"<<std::endl;
0064       return false;
0065     }
0066     matchingbranch=TdataIn->SetBranchAddress("event",&eventptr);
0067     if(matchingbranch<0){
0068       std::cout<<"Error retrieving Event info from the tree"<<std::endl;
0069       return false;
0070     }
0071     
0072     TcalibIn = (TTree*) RootInput->Get("Calib");
0073     if(!TcalibIn){
0074       std::cout<<"Could not retrieve Calib tree, leaving"<<std::endl;
0075       //return false;
0076     }
0077     else {
0078       matchingbranch=TcalibIn->SetBranchAddress("calib",&calibptr);
0079       if(matchingbranch<0){
0080         std::cout<<"Error retrieving calibration info from the tree"<<std::endl;
0081         TcalibIn=nullptr;
0082       }
0083     }
0084   } else {
0085     std::cout<<"An input file is required, aborting"<<std::endl;
0086     return false;
0087   }  
0088   
0089   //Setup Output files
0090   if(!RootOutputName.IsNull()){    
0091     
0092     bool sCOF = CreateOutputRootFile();
0093     if (!sCOF) return false;
0094     
0095     TsetupOut = new TTree("Setup","Setup");
0096     setup=Setup::GetInstance();
0097     TsetupOut->Branch("setup",&rsw);
0098     TdataOut = new TTree("Data","Data");
0099     TdataOut->Branch("event",&event);
0100     TcalibOut = new TTree("Calib","Calib");
0101     TcalibOut->Branch("calib",&calib);
0102   } 
0103   
0104   //Setup Output files
0105   if(!RootOutputNameHist.IsNull()){    
0106     std::cout<<"Creating output file for hists: "<< RootOutputNameHist.Data() <<std::endl;
0107     bool sCOF = CreateOutputRootFileHist();
0108     if (!sCOF) return false;
0109   }
0110   
0111   if ( RootOutputName.IsNull() && RootOutputNameHist.IsNull()){
0112     std::cout<<"Output file name is missing, please add"<<std::endl;
0113     return false;
0114   }  
0115   
0116   return true;
0117 }
0118 
0119 // ****************************************************************************
0120 // Primary process function to start all calibrations
0121 // ****************************************************************************
0122 bool DataAnalysis::Process(void){
0123   bool status;
0124   ROOT::EnableImplicitMT();
0125   
0126   if(RunQA){
0127     status=QAData();
0128   }  
0129   
0130   if(RunSimpleQA){
0131     status=SimpleQAData();
0132   }  
0133   return status;
0134 }
0135 
0136 
0137 
0138 //***********************************************************************************************
0139 //*********************** Advanced QA routine ***************************************************
0140 //***********************************************************************************************
0141 bool DataAnalysis::QAData(void){
0142   std::cout<<"QA data"<<std::endl;
0143   if(debug==1000){
0144     std::cerr<<"Time Min: "<<timemin<<std::endl;
0145     std::cerr<<"Time Max: "<<timemax<<std::endl;
0146   }
0147   std::map<int,RunInfo> ri=readRunInfosFromFile(RunListInputName.Data(),debug,0);
0148   TdataIn->GetEntry(0);
0149   Int_t runNr           = event.GetRunNumber();
0150   ReadOut::Type typeRO  = event.GetROtype();
0151 
0152   // Get run info object
0153   std::map<int,RunInfo>::iterator it=ri.find(runNr);
0154 
0155   int species = -1;
0156   species = GetSpeciesIntFromRunInfo(it->second);
0157   float beamE = it->second.energy;
0158   std::cout << "Beam energy:" << beamE << std::endl;
0159   if (species == -1){
0160       std::cout << "WARNING: species unknown: " << it->second.species.Data() << "  aborting"<< std::endl;
0161       return false;
0162   }
0163   
0164   Int_t maxCellsPerLayer = setup->GetMaxChannelInLayerFull()+1;
0165   //=============================================================================================
0166   // Create output histos
0167   //=============================================================================================
0168   // create HG and LG histo's per channel
0169   TH2D* hspectraHGCorrvsCellID      = nullptr;
0170   TH2D* hspectraLGCorrvsCellID      = nullptr;
0171   TH2D* hspectraHGvsCellID          = nullptr;
0172   TH2D* hspectraLGvsCellID          = nullptr;
0173   TH2D* hspectraTOTvsCellID         = nullptr;
0174 
0175   if (typeRO == ReadOut::Type::Caen) {
0176     hspectraHGCorrvsCellID      = new TH2D( "hspectraHGCorr_vsCellID","ADC spectrum High Gain corrected vs CellID; cell ID; ADC_{HG} (arb. units)  ; counts ",
0177                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5, 4000,-200,3800);
0178     hspectraHGCorrvsCellID->SetDirectory(0);
0179     hspectraLGCorrvsCellID      = new TH2D( "hspectraLGCorr_vsCellID","ADC spectrum Low Gain corrected vs CellID; cell ID; ADC_{LG} (arb. units)  ; counts  ",
0180                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5, 4000,-200,3800);
0181     hspectraLGCorrvsCellID->SetDirectory(0);
0182     hspectraHGvsCellID      = new TH2D( "hspectraHG_vsCellID","ADC spectrum High Gain vs CellID; cell ID; ADC_{HG} (arb. units)   ; counts ",
0183                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5, 4000,0,4000);
0184     hspectraHGvsCellID->SetDirectory(0);
0185     hspectraLGvsCellID      = new TH2D( "hspectraLG_vsCellID","ADC spectrum Low Gain vs CellID; cell ID; ADC_{LG} (arb. units)  ; counts",
0186                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5, 4000,0,4000);
0187     hspectraLGvsCellID->SetDirectory(0);
0188   } else if (typeRO == ReadOut::Type::Hgcroc) {
0189     hspectraHGvsCellID      = new TH2D( "hspectraADCvsCellID","ADC vs CellID; cell ID; ADC (arb. units)  ; counts ",
0190                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5, 1128,-100,1028);
0191     hspectraHGvsCellID->SetDirectory(0);
0192     hspectraTOTvsCellID= new TH2D( "hspectraTOTvsCellID","TOT vs CellID; cell ID; TOT (arb. units)  ; counts ",
0193                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5, 4096,0,4096);
0194     hspectraTOTvsCellID->SetDirectory(0);
0195   }
0196   TH2D* hspectraEnergyvsCellID  = new TH2D( "hspectraEnergy_vsCellID","Energy vs CellID; cell ID; E (mip eq./tile)    ; counts",
0197                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5, 8000,0,250);
0198   hspectraEnergyvsCellID->SetDirectory(0);
0199 
0200   // beam species dependent binning for energy hists:
0201   TH2D* hspectraEnergyTotvsNCells     = nullptr;
0202   TH2D* hspectraEnergyTotvsNCellsMuon = nullptr;
0203   TH2D* hspectraEnergyVsLayer         = nullptr;
0204   TH2D* hspectraEnergyVsLayerMuon     = nullptr;
0205   // muons 
0206   if (species == 0 || species == 2 ){
0207     hspectraEnergyTotvsNCells       = new TH2D( "hspectraTotEnergy_vsNCells","Energy vs CellID; N_{Cells}; E_{tot} (mip eq./tile) ; counts",
0208                                                 setup->GetNActiveCells()+1, -0.5, setup->GetNActiveCells()+1-0.5, 4000,0,500);
0209     hspectraEnergyTotvsNCellsMuon   = new TH2D( "hspectraTotEnergy_vsNCellsMuon","Energy vs CellID; N_{Cells}; E_{tot} (mip eq./tile) ; counts",
0210                                                 setup->GetNActiveCells()+1, -0.5, setup->GetNActiveCells()+1-0.5, 4000,0,500);    
0211     hspectraEnergyVsLayer           = new TH2D( "hspectraTotLayerEnergy_vsLayer","Energy in layer vs Layer; Layer; E_{layer} (mip eq./tile) ; counts",
0212                                                 setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, 2000,0,400);
0213     hspectraEnergyVsLayerMuon       = new TH2D( "hspectraTotLayerEnergy_vsLayerMuon","Energy in layer vs Layer Muon triggers; Layer; E_{layer} (mip eq./tile) ; counts",
0214                                                 setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, 2000,0,400);
0215   // em-probes
0216   } else if ( species == 1 || species == 3 || species == 6 ){
0217     hspectraEnergyTotvsNCells       = new TH2D( "hspectraTotEnergy_vsNCells","Energy vs CellID; N_{Cells}; E_{tot} (mip eq./tile) ; counts",
0218                                                 setup->GetNActiveCells()+1, -0.5, setup->GetNActiveCells()+1-0.5, 8000,0,1000);
0219     hspectraEnergyTotvsNCellsMuon   = new TH2D( "hspectraTotEnergy_vsNCellsMuon","Energy vs CellID; N_{Cells}; E_{tot} (mip eq./tile) ; counts",
0220                                                 setup->GetNActiveCells()+1, -0.5, setup->GetNActiveCells()+1-0.5, 8000,0,1000);    
0221     hspectraEnergyVsLayer           = new TH2D( "hspectraTotLayerEnergy_vsLayer","Energy in layer vs Layer; Layer; E_{layer} (mip eq./tile) ; counts",
0222                                                 setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, 2000,0,400);
0223     hspectraEnergyVsLayerMuon       = new TH2D( "hspectraTotLayerEnergy_vsLayerMuon","Energy in layer vs Layer Muon triggers; Layer; E_{layer} (mip eq./tile) ; counts",
0224                                                 setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, 2000,0,400);
0225   // hadrons
0226   } else if ( species == 4 || species == 5  ){
0227     hspectraEnergyTotvsNCells       = new TH2D( "hspectraTotEnergy_vsNCells","Energy vs CellID; N_{Cells}; E_{tot} (mip eq./tile) ; counts",
0228                                                 setup->GetNActiveCells()+1, -0.5, setup->GetNActiveCells()+1-0.5, 16000,0,2000);
0229     hspectraEnergyTotvsNCellsMuon   = new TH2D( "hspectraTotEnergy_vsNCellsMuon","Energy vs CellID; N_{Cells}; E_{tot} (mip eq./tile) ; counts",
0230                                                 setup->GetNActiveCells()+1, -0.5, setup->GetNActiveCells()+1-0.5, 16000,0,2000);
0231     hspectraEnergyVsLayer           = new TH2D( "hspectraTotLayerEnergy_vsLayer","Energy in layer vs Layer; Layer; E_{layer} (mip eq./tile) ; counts",
0232                                                 setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, 4000,0,800);
0233     hspectraEnergyVsLayerMuon       = new TH2D( "hspectraTotLayerEnergy_vsLayerMuon","Energy in layer vs Layer Muon triggers; Layer; E_{layer} (mip eq./tile) ; counts",
0234                                                 setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, 4000,0,800);
0235   }
0236   
0237   hspectraEnergyTotvsNCells->SetDirectory(0);
0238   hspectraEnergyTotvsNCells->Sumw2();
0239   hspectraEnergyTotvsNCellsMuon->SetDirectory(0);
0240   hspectraEnergyTotvsNCellsMuon->Sumw2();
0241   hspectraEnergyVsLayer->SetDirectory(0);
0242   hspectraEnergyVsLayer->Sumw2();
0243   hspectraEnergyVsLayerMuon->SetDirectory(0);
0244   hspectraEnergyVsLayerMuon->Sumw2();
0245   
0246   Double_t minXPos = setup->GetMinX();
0247   Double_t maxXPos = setup->GetMaxX();
0248   Double_t minYPos = setup->GetMinY();
0249   Double_t maxYPos = setup->GetMaxY();
0250   Int_t nBinsXPos  = (Int_t)((maxXPos-minXPos)/0.1);
0251   Int_t nBinsYPos  = (Int_t)((maxXPos-minXPos)/0.1);
0252   
0253   TH2D* hAverageXVsLayer  = new TH2D( "hAverageX_vsLayer","Av. X pos in layer vs Layer; Layer; X_{pos} (cm) ; counts",
0254                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, nBinsXPos,minXPos,maxXPos);
0255   hAverageXVsLayer->SetDirectory(0);
0256   hAverageXVsLayer->Sumw2();
0257   TH2D* hAverageXVsLayerMuon  = new TH2D( "hAverageX_vsLayerMuon","Av. X pos in layer vs Layer Muon triggers; Layer; X_{pos} (cm) ; counts",
0258                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, nBinsXPos,minXPos,maxXPos);
0259   hAverageXVsLayerMuon->SetDirectory(0);
0260   hAverageXVsLayer->Sumw2();
0261   TH2D* hAverageYVsLayer  = new TH2D( "hAverageX_vsLayer","Av. Y pos in layer vs Layer; Layer; Y_{pos} (cm) ; counts",
0262                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, nBinsYPos,minYPos,maxYPos);
0263   hAverageYVsLayer->SetDirectory(0);
0264   hAverageYVsLayer->Sumw2();
0265   TH2D* hAverageYVsLayerMuon  = new TH2D( "hAverageX_vsLayerMuon","Av. Y pos in layer vs Layer Muon triggers; Layer; Y_{pos} (cm) ; counts",
0266                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, nBinsYPos,minYPos,maxYPos);
0267   hAverageYVsLayerMuon->SetDirectory(0);
0268   hAverageYVsLayerMuon->Sumw2();
0269   TH2D* hNCellsVsLayer  = new TH2D( "hnCells_vsLayer","NCells in layer vs Layer; Layer; N_{cells,layer} ; counts",
0270                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxCellsPerLayer,-0.5,maxCellsPerLayer-0.5);
0271   hNCellsVsLayer->SetDirectory(0);
0272   hNCellsVsLayer->Sumw2();
0273   TH2D* hNCellsVsLayerMuon  = new TH2D( "hnCells_vsLayerMuon","NCells in layer vs Layer  Muon triggers; Layer; N_{cells,layer} ; counts",
0274                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxCellsPerLayer,-0.5,maxCellsPerLayer-0.5);
0275   hNCellsVsLayerMuon->SetDirectory(0);
0276   hNCellsVsLayerMuon->Sumw2();
0277   
0278 
0279 
0280   //Create 1D Histos for Delta time
0281   TH1D* hDeltaTime = new TH1D("hDeltaTime", "Time Difference between Events; Delta Time (#mus); Counts", 2000, 0, 100000);
0282   
0283   std::map<int,TileSpectra> hSpectra;
0284   std::map<int,TileSpectra> hSpectraTrigg;
0285   std::map<int, TileSpectra>::iterator ithSpectra;
0286   std::map<int, TileSpectra>::iterator ithSpectraTrigg;
0287 
0288   RootOutputHist->mkdir("IndividualCells");
0289   RootOutputHist->mkdir("IndividualCellsTrigg");
0290   
0291   //=============================================================================================
0292   // Setup general variables for analysis
0293   //=============================================================================================
0294   std::cout << "starting to run QA"<< std::endl;
0295   // initialize calib
0296   TcalibIn->GetEntry(0);
0297   int actCh1st = 0;
0298   double averageScale = calib.GetAverageScaleHigh(actCh1st);
0299   std::cout << "average HG mip: " << averageScale << "\t active ch: "<< actCh1st<< std::endl;
0300 
0301   // get and set maximum event number
0302   int evts=TdataIn->GetEntries();
0303   if ((eventNumber > -1) && (eventNumber <= evts)){
0304     evts = eventNumber;
0305     std::cout<<"restricting number of events in QA to " << evts << std::endl;
0306   }
0307   int evtsMuon= 0;
0308   double last_time = 0;
0309   double DeltaTime = 1000000000;
0310   double minESave  = 0.;
0311   
0312   int nActiveLayers = setup->GetNActiveLayers();
0313   
0314   //=============================================================================================
0315   // Run over all events in tree
0316   //=============================================================================================  
0317   for(int i=0; i<evts; i++){
0318     TdataIn->GetEntry(i);
0319     if (i%5000 == 0 && debug > 0) std::cout << "Reading " <<  i << " / " << evts << " events" << std::endl;
0320     
0321     //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
0322     // Set up event variables
0323     //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
0324     double Etot = 0;
0325     int nCells  = 0;
0326     bool muontrigg = false;
0327     int locMuon = 0;
0328     std::map<int,Layer> layers;
0329     std::map<int, Layer>::iterator ithLayer;
0330     
0331     // Find and fill delta time
0332     double current_time = event.GetTimeStamp();
0333     if(last_time != 0){
0334       DeltaTime = current_time - last_time;
0335       if (debug == 1000){
0336         std::cerr<< "current timestamp: "<< current_time<<std::endl;
0337       }
0338     }
0339     if (DeltaTime == 0){
0340       if(debug == 1001){
0341         std::cerr<< "Run Number: " << runNr <<"      Event Number: "<< i << std::endl;
0342         std::cerr<< "Previous Timestamp (us): "<< last_time <<"      Current Timestamp: "<< current_time<<std::endl;
0343       }
0344     }
0345         //cut based on DeltaTime range
0346     if(DeltaTime < timemin || DeltaTime > timemax){
0347       //std::cout<<"event rejected:"<< i <<std::endl;
0348       last_time = current_time;
0349       continue;
0350     }
0351     hDeltaTime->Fill(DeltaTime);
0352     last_time = current_time;
0353         
0354     //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
0355     // Evaluate all cells in event
0356     //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
0357     for(int j=0; j<event.GetNTiles(); j++){
0358       if (typeRO == ReadOut::Type::Caen) {
0359         Caen* aTile=(Caen*)event.GetTile(j);
0360         double corrHG = aTile->GetADCHigh()-calib.GetPedestalMeanH(aTile->GetCellID());
0361         double corrLG = aTile->GetADCLow()-calib.GetPedestalMeanL(aTile->GetCellID());
0362         hspectraHGvsCellID->Fill(aTile->GetCellID(), aTile->GetADCHigh());
0363         hspectraLGvsCellID->Fill(aTile->GetCellID(), aTile->GetADCLow());
0364         hspectraHGCorrvsCellID->Fill(aTile->GetCellID(), corrHG);
0365         hspectraLGCorrvsCellID->Fill(aTile->GetCellID(), corrLG);
0366         if(aTile->GetE()> minESave ){ 
0367           nCells++;
0368           int currLayer = setup->GetLayer(aTile->GetCellID());
0369           ithLayer=layers.find(currLayer);
0370           if(ithLayer!=layers.end()){
0371             ithLayer->second.nCells+=1;
0372             ithLayer->second.energy+=aTile->GetE();
0373             ithLayer->second.avX+=setup->GetX(aTile->GetCellID())*aTile->GetE();
0374             ithLayer->second.avY+=setup->GetY(aTile->GetCellID())*aTile->GetE();
0375           } else {
0376             layers[currLayer]=Layer();
0377             layers[currLayer].nCells+=1;
0378             layers[currLayer].energy+=aTile->GetE();
0379             layers[currLayer].avX+=setup->GetX(aTile->GetCellID())*aTile->GetE();
0380             layers[currLayer].avY+=setup->GetY(aTile->GetCellID())*aTile->GetE();
0381           }
0382         }
0383         // read current tile
0384         ithSpectraTrigg=hSpectraTrigg.find(aTile->GetCellID());
0385         if (aTile->GetLocalTriggerBit() == 1){
0386           if(ithSpectraTrigg!=hSpectraTrigg.end()){
0387             ithSpectraTrigg->second.FillTrigger(aTile->GetLocalTriggerPrimitive());
0388             ithSpectraTrigg->second.FillSpectraCAEN(corrLG,corrHG);
0389           } else {
0390             RootOutputHist->cd("IndividualCellsTrigg");
0391             hSpectraTrigg[aTile->GetCellID()]=TileSpectra("muonTriggers",aTile->GetCellID(),calib.GetTileCalib(aTile->GetCellID()),event.GetROtype(),debug);
0392             hSpectraTrigg[aTile->GetCellID()].FillTrigger(aTile->GetLocalTriggerPrimitive());
0393             hSpectraTrigg[aTile->GetCellID()].FillSpectraCAEN(corrLG,corrHG);
0394           }
0395           locMuon++;        
0396         }      
0397         
0398         ithSpectra=hSpectra.find(aTile->GetCellID());
0399         if (ithSpectra!=hSpectra.end()){
0400           ithSpectra->second.FillSpectraCAEN(corrLG,corrHG);
0401           ithSpectra->second.FillTrigger(aTile->GetLocalTriggerPrimitive());;
0402           ithSpectra->second.FillCorrCAEN(corrLG,corrHG);
0403         } else {
0404           RootOutputHist->cd("IndividualCells");
0405           hSpectra[aTile->GetCellID()]=TileSpectra("AllTriggers",aTile->GetCellID(),calib.GetTileCalib(aTile->GetCellID()),event.GetROtype(),debug);
0406           hSpectra[aTile->GetCellID()].FillSpectraCAEN(corrLG,corrHG);;
0407           hSpectra[aTile->GetCellID()].FillTrigger(aTile->GetLocalTriggerPrimitive());;
0408           hSpectra[aTile->GetCellID()].FillCorrCAEN(corrLG,corrHG);
0409         }
0410       } else if (typeRO == ReadOut::Type::Hgcroc){
0411         Hgcroc* aTile=(Hgcroc*)event.GetTile(j);
0412         double adc = aTile->GetIntegratedADC();
0413         double tot = aTile->GetRawTOT();
0414         double toa = aTile->GetRawTOA();
0415         hspectraHGvsCellID->Fill(aTile->GetCellID(), adc);
0416         hspectraTOTvsCellID->Fill(aTile->GetCellID(), tot);
0417         if(aTile->GetE() > minESave ){ 
0418           nCells++;
0419           int currLayer = setup->GetLayer(aTile->GetCellID());
0420           ithLayer=layers.find(currLayer);
0421           if(ithLayer!=layers.end()){
0422             ithLayer->second.nCells+=1;
0423             ithLayer->second.energy+=aTile->GetE();
0424             ithLayer->second.avX+=setup->GetX(aTile->GetCellID())*aTile->GetE();
0425             ithLayer->second.avY+=setup->GetY(aTile->GetCellID())*aTile->GetE();
0426           } else {
0427             layers[currLayer]=Layer();
0428             layers[currLayer].nCells+=1;
0429             layers[currLayer].energy+=aTile->GetE();
0430             layers[currLayer].avX+=setup->GetX(aTile->GetCellID())*aTile->GetE();
0431             layers[currLayer].avY+=setup->GetY(aTile->GetCellID())*aTile->GetE();
0432           }
0433         }
0434         // read current tile
0435         ithSpectraTrigg=hSpectraTrigg.find(aTile->GetCellID());
0436         if (aTile->GetLocalTriggerBit() == 1){
0437           if(ithSpectraTrigg!=hSpectraTrigg.end()){
0438             ithSpectraTrigg->second.FillTrigger(aTile->GetLocalTriggerPrimitive());
0439             ithSpectraTrigg->second.FillHGCROC(adc,toa,tot);
0440           } else {
0441             RootOutputHist->cd("IndividualCellsTrigg");
0442             hSpectraTrigg[aTile->GetCellID()]=TileSpectra("muonTriggers",aTile->GetCellID(),calib.GetTileCalib(aTile->GetCellID()),event.GetROtype(),debug);
0443             hSpectraTrigg[aTile->GetCellID()].FillTrigger(aTile->GetLocalTriggerPrimitive());
0444             hSpectraTrigg[aTile->GetCellID()].FillHGCROC(adc,toa,tot);
0445           }
0446           locMuon++;        
0447         }      
0448         
0449         ithSpectra=hSpectra.find(aTile->GetCellID());
0450         if (ithSpectra!=hSpectra.end()){
0451           ithSpectra->second.FillTrigger(aTile->GetLocalTriggerPrimitive());;
0452           ithSpectra->second.FillHGCROC(adc,toa,tot);
0453         } else {
0454           RootOutputHist->cd("IndividualCells");
0455           hSpectra[aTile->GetCellID()]=TileSpectra("AllTriggers",aTile->GetCellID(),calib.GetTileCalib(aTile->GetCellID()),event.GetROtype(),debug);
0456           hSpectra[aTile->GetCellID()].FillTrigger(aTile->GetLocalTriggerPrimitive());;
0457           hSpectra[aTile->GetCellID()].FillHGCROC(adc,toa,tot);
0458         }
0459       }
0460     }
0461     if (nCells > 0) {
0462 
0463       int nLayerSingleCell = 0;
0464       for(ithLayer=layers.begin(); ithLayer!=layers.end(); ++ithLayer){
0465         if (ithLayer->second.nCells == 1) 
0466             nLayerSingleCell++;
0467       }
0468       double fracLayer1Cell = (double)nLayerSingleCell/(int)layers.size();
0469       // if ( fracLayer1Cell > 0.8){
0470       //   muontrigg = true;
0471       //   evtsMuon++;        
0472       // }
0473       double fracLocMuon = (double)locMuon/nCells;
0474       if (fracLocMuon > 0.6){ 
0475         muontrigg = true;
0476         evtsMuon++;
0477       }
0478       
0479       if (muontrigg && debug > 3){
0480           std::cout << "Muon triggered:\t" <<  fracLayer1Cell << "\t" << fracLocMuon << "\t" << nActiveLayers << "\t" << nCells<< std::endl;
0481       }
0482       
0483       for(int l = 0; l < setup->GetNMaxLayer()+1;l++){
0484         ithLayer=layers.find(l);
0485         if (ithLayer!=layers.end()){          
0486           double avx     = ithLayer->second.avX/ithLayer->second.energy;
0487           double avy     = ithLayer->second.avY/ithLayer->second.energy;
0488           hspectraEnergyVsLayer->Fill(l,ithLayer->second.energy);
0489           if (muontrigg){
0490             hspectraEnergyVsLayerMuon->Fill(l,ithLayer->second.energy);
0491             hAverageXVsLayerMuon->Fill(l,avx);
0492             hAverageYVsLayerMuon->Fill(l,avy);
0493             hNCellsVsLayerMuon->Fill(l,ithLayer->second.nCells);
0494           }
0495           hAverageXVsLayer->Fill(l,avx);
0496           hAverageYVsLayer->Fill(l,avy);
0497           hNCellsVsLayer->Fill(l,ithLayer->second.nCells);
0498           if (ithLayer->second.nCells == 1) 
0499             nLayerSingleCell++;
0500         } else {
0501           hspectraEnergyVsLayer->Fill(l,0.);
0502           if (muontrigg){
0503             hspectraEnergyVsLayerMuon->Fill(l,0.);
0504             hAverageXVsLayerMuon->Fill(l,-20.);
0505             hAverageYVsLayerMuon->Fill(l,-20.);
0506             hNCellsVsLayerMuon->Fill(l,0);            
0507           }
0508           hAverageXVsLayer->Fill(l,-20);
0509           hAverageYVsLayer->Fill(l,-20);
0510           hNCellsVsLayer->Fill(l,0);
0511         }
0512       }
0513 
0514       for(int j=0; j<event.GetNTiles(); j++){
0515         Caen* aTile=(Caen*)event.GetTile(j);
0516         // remove bad channels from output
0517         double energy = aTile->GetE();
0518         if(energy!=0){ 
0519           hspectraEnergyvsCellID->Fill(aTile->GetCellID(), energy);
0520           Etot=Etot+energy;
0521         } 
0522       }
0523 
0524       hspectraEnergyTotvsNCells->Fill(nCells,Etot);
0525       if(muontrigg) hspectraEnergyTotvsNCellsMuon->Fill(nCells,Etot);
0526     }
0527   }
0528   
0529   if (IsCalibSaveToFile()){
0530     TString fileCalibPrint = RootOutputName;
0531     fileCalibPrint         = fileCalibPrint.ReplaceAll(".root","_calib.txt");
0532     calib.PrintCalibToFile(fileCalibPrint);
0533   }
0534   
0535   RootInput->Close();      
0536   
0537   RootOutputHist->cd();
0538 
0539   hDeltaTime->Write();
0540 
0541   hspectraHGvsCellID->Write();
0542   if (typeRO == ReadOut::Type::Caen){
0543     hspectraLGvsCellID->Write();
0544     hspectraHGCorrvsCellID->Write();
0545     hspectraLGCorrvsCellID->Write();    
0546   } else if (typeRO == ReadOut::Type::Hgcroc){
0547     hspectraTOTvsCellID->Write();
0548   }
0549   
0550   hspectraEnergyvsCellID->Write();
0551   hspectraEnergyTotvsNCells->Write();
0552   hspectraEnergyTotvsNCellsMuon->Write();
0553   
0554   TH2D* hspectraEnergyTotvsNCells_WoMuon = (TH2D*)hspectraEnergyTotvsNCells->Clone("hspectraTotEnergy_vsNCells_woMuon");
0555   // hspectraEnergyTotvsNCells_WoMuon->Sumw2();
0556   hspectraEnergyTotvsNCells_WoMuon->Add(hspectraEnergyTotvsNCellsMuon, -1);
0557   hspectraEnergyTotvsNCells_WoMuon->Write();
0558   
0559   TH1D* hspectraEnergyTot = (TH1D*)hspectraEnergyTotvsNCells->ProjectionY();
0560   hspectraEnergyTot->SetDirectory(0);
0561   TH1D* hspectraNCells = (TH1D*)hspectraEnergyTotvsNCells->ProjectionX();
0562   hspectraNCells->SetDirectory(0);
0563   hspectraEnergyTot->Write("hTotEnergy");
0564   hspectraNCells->Write("hNCells");
0565 
0566   TH1D* hspectraEnergyTotMuon = (TH1D*)hspectraEnergyTotvsNCellsMuon->ProjectionY();
0567   hspectraEnergyTotMuon->SetDirectory(0);
0568   TH1D* hspectraNCellsMuon = (TH1D*)hspectraEnergyTotvsNCellsMuon->ProjectionX();
0569   hspectraNCellsMuon->SetDirectory(0);
0570   hspectraEnergyTotMuon->Write("hTotEnergyMuon");
0571   hspectraNCellsMuon->Write("hNCellsMuon");
0572   
0573   TH1D* hspectraEnergyTotNonMuon = (TH1D*)hspectraEnergyTotvsNCells_WoMuon->ProjectionY();
0574   hspectraEnergyTotNonMuon->SetDirectory(0);
0575   TH1D* hspectraNCellsNonMuon = (TH1D*)hspectraEnergyTotvsNCells_WoMuon->ProjectionX();
0576   hspectraNCellsNonMuon->SetDirectory(0);
0577   hspectraEnergyTotNonMuon->Write("hTotEnergyNonMuon");
0578   hspectraNCellsNonMuon->Write("hNCellsNonMuon");
0579   
0580   // position hists
0581   hAverageXVsLayer->Write();
0582   hAverageXVsLayerMuon->Write();
0583   TH2D* hAverageXVsLayer_WoMuon = (TH2D*)hAverageXVsLayer->Clone("hAverageX_vsLayer_woMuon");
0584   // hAverageXVsLayer_WoMuon->Sumw2();
0585   hAverageXVsLayer_WoMuon->Add(hAverageXVsLayerMuon, -1);
0586   hAverageXVsLayer_WoMuon->Write();
0587   
0588   hAverageYVsLayer->Write();
0589   hAverageYVsLayerMuon->Write();
0590   TH2D* hAverageYVsLayer_WoMuon = (TH2D*)hAverageYVsLayer->Clone("hAverageY_vsLayer_woMuon");
0591   // hAverageYVsLayer_WoMuon->Sumw2();
0592   hAverageYVsLayer_WoMuon->Add(hAverageYVsLayerMuon, -1);
0593   hAverageYVsLayer_WoMuon->Write();
0594   
0595   hNCellsVsLayer->Write();
0596   hNCellsVsLayerMuon->Write();
0597   TH2D* hNCellsVsLayer_WoMuon = (TH2D*)hNCellsVsLayer->Clone("hnCells_vsLayer_woMuon");
0598   // hNCellsVsLayer_WoMuon->Sumw2();
0599   hNCellsVsLayer_WoMuon->Add(hNCellsVsLayerMuon, -1);
0600   hNCellsVsLayer_WoMuon->Write();
0601   
0602   hspectraEnergyVsLayer->Write();
0603   // hspectraEnergyVsLayer->Sumw2();
0604   hspectraEnergyVsLayerMuon->Write();
0605   // hspectraEnergyVsLayerMuon->Sumw2();
0606   TH2D* hspectraEnergyVsLayer_WoMuon = (TH2D*)hspectraEnergyVsLayer->Clone("hspectraTotLayerEnergy_vsLayer_woMuon");
0607   // hspectraEnergyVsLayer_WoMuon->Sumw2();
0608   hspectraEnergyVsLayer_WoMuon->Add(hspectraEnergyVsLayerMuon, -1);
0609   hspectraEnergyVsLayer_WoMuon->Write();
0610 
0611   TH1D* histELayer_All[65];
0612   TH1D* histELayer_Muon[65];
0613   TH1D* histELayer_WoMuon[65];
0614   TH1D* histXLayer_All[65];
0615   TH1D* histXLayer_Muon[65];
0616   TH1D* histXLayer_WoMuon[65];
0617   TH1D* histYLayer_All[65];
0618   TH1D* histYLayer_Muon[65];
0619   TH1D* histYLayer_WoMuon[65];
0620   TH1D* histNCellsLayer_All[65];
0621   TH1D* histNCellsLayer_Muon[65];
0622   TH1D* histNCellsLayer_WoMuon[65];
0623   Float_t maxYLAll    = 0;
0624   Float_t maxYLMuon   = 0;
0625   Float_t maxYLWOMuon = 0;
0626   Float_t maxXLAll    = 0;
0627   Float_t maxXLMuon   = 0;
0628   Float_t maxXLWOMuon = 0;
0629   std::map<int,Layer> layersMeanAll;
0630   std::map<int, Layer>::iterator ithLayerMeanAll;
0631   std::map<int,Layer> layersMeanMuon;
0632   std::map<int, Layer>::iterator ithLayerMeanMuon;
0633   std::map<int,Layer> layersMeanWOMuon;
0634   std::map<int, Layer>::iterator ithLayerMeanWOMuon;
0635   
0636   TGraphErrors* graphMeanE_Layer        = new TGraphErrors();
0637   TGraphErrors* graphMeanE_Layer_muon   = new TGraphErrors();
0638   TGraphErrors* graphMeanE_Layer_woMuon = new TGraphErrors();
0639   
0640   for(int l = 0; l < setup->GetNMaxLayer()+1;l++){  
0641     
0642     bool layerOn = true;
0643     
0644     histNCellsLayer_All[l]     = (TH1D*)hNCellsVsLayer->ProjectionY(Form("histNCellsLayer_AllTrigg_Layer_%02d",l),l+1,l+1);
0645     histNCellsLayer_All[l]->Write();
0646     histNCellsLayer_Muon[l]    = (TH1D*)hNCellsVsLayerMuon->ProjectionY(Form("histNCellsLayer_MuonTrigg_Layer_%02d",l),l+1,l+1);
0647     histNCellsLayer_Muon[l]->Write();
0648     histNCellsLayer_WoMuon[l]  = (TH1D*)hNCellsVsLayer_WoMuon->ProjectionY(Form("histNCellsLayer_WOMuonTrigg_Layer_%02d",l),l+1,l+1);
0649     histNCellsLayer_WoMuon[l]->Write();
0650 
0651     
0652     if (histNCellsLayer_All[l]->GetMean() < 0.001) layerOn =false;
0653     
0654     histELayer_All[l]     = (TH1D*)hspectraEnergyVsLayer->ProjectionY(Form("histELayer_AllTrigg_Layer_%02d",l),l+1,l+1);
0655     // histELayer_All[l]->Sumw2();
0656     histELayer_All[l]->Write();
0657     histELayer_Muon[l]    = (TH1D*)hspectraEnergyVsLayerMuon->ProjectionY(Form("histELayer_MuonTrigg_Layer_%02d",l),l+1,l+1);
0658     // histELayer_Muon[l]->Sumw2();
0659     histELayer_Muon[l]->Write();
0660     histELayer_WoMuon[l]  = (TH1D*)hspectraEnergyVsLayer_WoMuon->ProjectionY(Form("histELayer_WOMuonTrigg_Layer_%02d",l),l+1,l+1);
0661     // histELayer_WoMuon[l]->Sumw2();
0662     histELayer_WoMuon[l]->Write();
0663     
0664     if (layerOn){
0665       if (maxYLAll < histELayer_All[l]->GetMaximum())maxYLAll =  histELayer_All[l]->GetMaximum();
0666       if (maxXLAll < FindLastBinXAboveMin(histELayer_All[l],2)) maxXLAll = FindLastBinXAboveMin(histELayer_All[l],2);
0667       if (maxYLMuon < histELayer_Muon[l]->GetMaximum() )maxYLMuon =  histELayer_Muon[l]->GetMaximum();
0668       if (maxXLMuon < FindLastBinXAboveMin(histELayer_Muon[l],2) ) maxXLMuon = FindLastBinXAboveMin(histELayer_Muon[l],2);
0669       if (maxYLWOMuon < histELayer_WoMuon[l]->GetMaximum() )maxYLWOMuon =  histELayer_WoMuon[l]->GetMaximum();
0670       if (maxXLWOMuon < FindLastBinXAboveMin(histELayer_WoMuon[l],2) ) maxXLWOMuon = FindLastBinXAboveMin(histELayer_WoMuon[l],2);
0671     }
0672 
0673     histXLayer_All[l]     = (TH1D*)hAverageXVsLayer->ProjectionY(Form("histXLayer_AllTrigg_Layer_%02d",l),l+1,l+1);
0674     histXLayer_All[l]->Write();
0675     histXLayer_Muon[l]    = (TH1D*)hAverageXVsLayerMuon->ProjectionY(Form("histXLayer_MuonTrigg_Layer_%02d",l),l+1,l+1);
0676     histXLayer_Muon[l]->Write();
0677     histXLayer_WoMuon[l]  = (TH1D*)hAverageXVsLayer_WoMuon->ProjectionY(Form("histXLayer_WOMuonTrigg_Layer_%02d",l),l+1,l+1);
0678     histXLayer_WoMuon[l]->Write();
0679 
0680     histYLayer_All[l]     = (TH1D*)hAverageYVsLayer->ProjectionY(Form("histYLayer_AllTrigg_Layer_%02d",l),l+1,l+1);
0681     histYLayer_All[l]->Write();
0682     histYLayer_Muon[l]    = (TH1D*)hAverageYVsLayerMuon->ProjectionY(Form("histYLayer_MuonTrigg_Layer_%02d",l),l+1,l+1);
0683     histYLayer_Muon[l]->Write();
0684     histYLayer_WoMuon[l]  = (TH1D*)hAverageYVsLayer_WoMuon->ProjectionY(Form("histYLayer_WOMuonTrigg_Layer_%02d",l),l+1,l+1);
0685     histYLayer_WoMuon[l]->Write();
0686 
0687     layersMeanAll[l]=Layer();
0688     layersMeanAll[l].energy+=histELayer_All[l]->GetMean();
0689     layersMeanAll[l].avX+=histXLayer_All[l]->GetMean();
0690     layersMeanAll[l].avY+=histYLayer_All[l]->GetMean();
0691     layersMeanAll[l].nCells+=histNCellsLayer_All[l]->GetMean();
0692     
0693     layersMeanMuon[l]=Layer();
0694     layersMeanMuon[l].energy+=histELayer_Muon[l]->GetMean();
0695     layersMeanMuon[l].avX+=histXLayer_Muon[l]->GetMean();
0696     layersMeanMuon[l].avY+=histYLayer_Muon[l]->GetMean();
0697     layersMeanMuon[l].nCells+=histNCellsLayer_Muon[l]->GetMean();
0698     
0699     layersMeanWOMuon[l]=Layer();
0700     layersMeanWOMuon[l].energy+=histELayer_WoMuon[l]->GetMean();
0701     layersMeanWOMuon[l].avX+=histXLayer_WoMuon[l]->GetMean();
0702     layersMeanWOMuon[l].avY+=histYLayer_WoMuon[l]->GetMean();
0703     layersMeanWOMuon[l].nCells+=histNCellsLayer_WoMuon[l]->GetMean();
0704     
0705     graphMeanE_Layer->SetPoint(graphMeanE_Layer->GetN(), l,histELayer_All[l]->GetMean());
0706     graphMeanE_Layer->SetPointError(graphMeanE_Layer->GetN()-1, 0,histELayer_All[l]->GetMeanError());
0707     // graphMeanE_Layer->SetPointError(graphMeanE_Layer->GetN()-1, 0,histELayer_All[l]->GetRMS());
0708     graphMeanE_Layer_muon->SetPoint(graphMeanE_Layer_muon->GetN(),l,histELayer_Muon[l]->GetMean());
0709     graphMeanE_Layer_muon->SetPointError(graphMeanE_Layer_muon->GetN()-1, 0,histELayer_Muon[l]->GetMeanError());
0710     // graphMeanE_Layer_muon->SetPointError(graphMeanE_Layer_muon->GetN()-1, 0,histELayer_Muon[l]->GetRMS());
0711     graphMeanE_Layer_woMuon->SetPoint(graphMeanE_Layer_woMuon->GetN(),l,histELayer_WoMuon[l]->GetMean());
0712     graphMeanE_Layer_woMuon->SetPointError(graphMeanE_Layer_woMuon->GetN()-1, 0,histELayer_WoMuon[l]->GetMeanError());
0713     // graphMeanE_Layer_woMuon->SetPointError(graphMeanE_Layer_woMuon->GetN()-1, 0,histELayer_WoMuon[l]->GetRMS());
0714   }
0715   
0716   std::cout<< "average per layer - All" << std::endl;
0717   graphMeanE_Layer->Write("graphMeanE_Layer");
0718   std::cout<< "average per layer - no Muon" << std::endl;
0719   graphMeanE_Layer_woMuon->Write("graphMeanE_Layer_woMuon");
0720   std::cout<< "average per layer - Muon" << std::endl;
0721   graphMeanE_Layer_muon->Write("graphMeanE_Layer_Muon");
0722   
0723   //**********************************************************************
0724   //********************* Plotting ***************************************
0725   //**********************************************************************  
0726   TString outputDirPlots = GetPlotOutputDir();
0727   gSystem->Exec("mkdir -p "+outputDirPlots);
0728   
0729   
0730   
0731   //**********************************************************************
0732   // Create canvases for channel overview plotting
0733   //**********************************************************************
0734   Double_t textSizeRel = 0.035;
0735   StyleSettingsBasics("pdf");
0736   SetPlotStyle();
0737   
0738   TCanvas* canvas2DCorr = new TCanvas("canvasCorrPlots","",0,0,1450,1300);  // gives the page size
0739   DefaultCanvasSettings( canvas2DCorr, 0.08, 0.13, 0.045, 0.07);
0740   canvas2DCorr->SetLogz(1);
0741   TCanvas* canvas2DCorrWOLine = new TCanvas("canvasCorrPlotsWoLine","",0,0,1450,1300);  // gives the page size
0742   DefaultCanvasSettings( canvas2DCorrWOLine, 0.08, 0.13, 0.01, 0.07);
0743   canvas2DCorrWOLine->SetLogz(1);
0744   
0745   TCanvas* canvasDeltaTime = new TCanvas("canvasDeltaTime","",0,0,1450,1300);  // gives the page size
0746   DefaultCanvasSettings( canvasDeltaTime, 0.08, 0.07, 0.01, 0.07);
0747   canvasDeltaTime->SetLogy(1);
0748   if(DeltaTimePlot>0){
0749     PlotSimple1D( canvasDeltaTime, hDeltaTime, -10000, timemax, textSizeRel, Form("%s/deltaTime.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1);
0750   } 
0751 
0752   PlotSimple2D( canvas2DCorr, hspectraHGvsCellID, -10000, -10000, textSizeRel, Form("%s/HG.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
0753   if (typeRO == ReadOut::Type::Caen){
0754     PlotSimple2D( canvas2DCorr, hspectraLGvsCellID, -10000, -10000, textSizeRel, Form("%s/LG.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
0755     PlotSimple2D( canvas2DCorr, hspectraHGCorrvsCellID, -10000, -10000, textSizeRel, Form("%s/HGCorr.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
0756     PlotSimple2D( canvas2DCorr, hspectraLGCorrvsCellID, -10000, -10000, textSizeRel, Form("%s/LGCorr.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
0757   } else if (typeRO == ReadOut::Type::Hgcroc){
0758     PlotSimple2D( canvas2DCorr, hspectraTOTvsCellID, -10000, -10000, textSizeRel, Form("%s/TOT.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);    
0759   }
0760   PlotSimple2D( canvas2DCorr, hspectraEnergyvsCellID, -10000, -10000, textSizeRel, Form("%s/EnergyVsCellID.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
0761   PlotSimple2D( canvas2DCorr, hspectraEnergyTotvsNCells, -10000, -10000, textSizeRel, Form("%s/EnergyTotalVsNCells.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true, Form("evts = %d",evts));
0762   PlotSimple2D( canvas2DCorr, hspectraEnergyTotvsNCellsMuon, -10000, -10000, textSizeRel, Form("%s/EnergyTotalVsNCells_MuonTrigg.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true, Form("Muon evts = %d",evtsMuon));
0763   PlotSimple2D( canvas2DCorr, hspectraEnergyTotvsNCells_WoMuon, -10000, -10000, textSizeRel, Form("%s/EnergyTotalVsNCells_WoMuonTrigg.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true, Form("Non Muon evts = %d",evts-evtsMuon));
0764 
0765   PlotSimple2D( canvas2DCorr, hspectraEnergyVsLayer, -10000, -10000, textSizeRel, Form("%s/EnergyVsLayer.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
0766   PlotSimple2D( canvas2DCorr, hspectraEnergyVsLayerMuon, -10000, -10000, textSizeRel, Form("%s/EnergyVsLayer_MuonTrigg.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true, "Muon triggers");
0767   PlotSimple2D( canvas2DCorr, hspectraEnergyVsLayer_WoMuon, -10000, -10000, textSizeRel, Form("%s/EnergyVsLayer_WOMuonTrigg.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true, "Non Muon triggers");
0768   PlotSimple2D( canvas2DCorr, hAverageXVsLayer, -10000, -10000, textSizeRel, Form("%s/XPosVsLayer.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
0769   PlotSimple2D( canvas2DCorr, hAverageXVsLayerMuon, -10000, -10000, textSizeRel, Form("%s/XPosVsLayer_MuonTrigg.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true, "Muon triggers");
0770   PlotSimple2D( canvas2DCorr, hAverageXVsLayer_WoMuon, -10000, -10000, textSizeRel, Form("%s/XPosVsLayer_WOMuonTrigg.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true, "Non Muon triggers");
0771   PlotSimple2D( canvas2DCorr, hAverageYVsLayer, -10000, -10000, textSizeRel, Form("%s/YPosVsLayer.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
0772   PlotSimple2D( canvas2DCorr, hAverageYVsLayerMuon, -10000, -10000, textSizeRel, Form("%s/YPosVsLayer_MuonTrigg.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true, "Muon triggers");
0773   PlotSimple2D( canvas2DCorr, hAverageYVsLayer_WoMuon, -10000, -10000, textSizeRel, Form("%s/YPosVsLayer_WOMuonTrigg.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true, "Non Muon triggers");
0774   PlotSimple2D( canvas2DCorr, hNCellsVsLayer, -10000, -10000, textSizeRel, Form("%s/NcellsLayerVsLayer.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
0775   PlotSimple2D( canvas2DCorr, hNCellsVsLayerMuon, -10000, -10000, textSizeRel, Form("%s/NcellsLayerVsLayer_MuonTrigg.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true, "Muon triggers");
0776   PlotSimple2D( canvas2DCorr, hNCellsVsLayer_WoMuon, -10000, -10000, textSizeRel, Form("%s/NcellsLayerVsLayer_WOMuonTrigg.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true, "Non Muon triggers");
0777   
0778   double maxEProp = 400.;
0779   if (TMath::Abs(beamE -5.0) < 0.01) {
0780       maxEProp = 300;
0781       std::cout << "reset max for beam E" << beamE << "\t" << maxEProp << std::endl;
0782   }
0783   if (species == 0 || species == 2 ){
0784     Plot2DWithGraph( canvas2DCorrWOLine, hspectraEnergyVsLayerMuon, nullptr, maxEProp, -10000,  textSizeRel, Form("%s/EnergyVsLayer_Propaganda.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", false, "");
0785     Plot2DWithGraph( canvas2DCorrWOLine, hspectraEnergyVsLayerMuon, graphMeanE_Layer_muon, maxEProp, -10000,  textSizeRel, Form("%s/EnergyVsLayer_PropagandaWGraph.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", false, "");
0786   } else {
0787     Plot2DWithGraph( canvas2DCorrWOLine, hspectraEnergyVsLayer_WoMuon,  nullptr, maxEProp, -10000, textSizeRel, Form("%s/EnergyVsLayer_Propaganda.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", false, "");
0788     Plot2DWithGraph( canvas2DCorrWOLine, hspectraEnergyVsLayer_WoMuon,  graphMeanE_Layer_woMuon, maxEProp, -10000, textSizeRel, Form("%s/EnergyVsLayer_PropagandaWGraph.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", false, "");
0789   }
0790   canvas2DCorrWOLine->SetLogy();
0791   if (species == 0 || species == 2 ){
0792     Plot2DWithGraph( canvas2DCorrWOLine, hspectraEnergyVsLayerMuon, nullptr, 400, -10000,  textSizeRel, Form("%s/EnergyVsLayer_Propaganda_LogY.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", false, "");
0793     Plot2DWithGraph( canvas2DCorrWOLine, hspectraEnergyVsLayerMuon, graphMeanE_Layer_muon, 400, -10000,  textSizeRel, Form("%s/EnergyVsLayer_PropagandaWGraph_LogY.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", false, "");
0794   } else {
0795     Plot2DWithGraph( canvas2DCorrWOLine, hspectraEnergyVsLayer_WoMuon,  nullptr, 400, -10000, textSizeRel, Form("%s/EnergyVsLayer_Propaganda_LogY.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", false, "");
0796     Plot2DWithGraph( canvas2DCorrWOLine, hspectraEnergyVsLayer_WoMuon,  graphMeanE_Layer_woMuon, 400, -10000, textSizeRel, Form("%s/EnergyVsLayer_PropagandaWGraph_LogY.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", false, "");
0797   }
0798   canvas2DCorrWOLine->SetLogy(kFALSE);
0799   
0800 
0801   TCanvas* canvas1DSimple = new TCanvas("canvas1DSimple","",0,0,1450,1300);  // gives the page size
0802   DefaultCanvasSettings( canvas1DSimple, 0.08, 0.01, 0.03, 0.07);
0803   hspectraEnergyTot->Rebin(4);
0804   hspectraEnergyTot->Scale(1./evts);
0805   hspectraEnergyTot->GetYaxis()->SetTitle("counts/event");
0806   
0807   TF1* fitGaussETot = new TF1("fitEtot","gaus",hspectraEnergyTot->GetMean()-1.5*hspectraEnergyTot->GetRMS(),hspectraEnergyTot->GetMean()+1.5*hspectraEnergyTot->GetRMS());
0808   hspectraEnergyTot->Fit(fitGaussETot,"QRMEN0S");
0809   std::cout<< "All trigg: "<<fitGaussETot->GetParameter(1) << "\t" << fitGaussETot->GetParameter(2) << std::endl;
0810   
0811   Double_t maxEPlot =  FindLastBinXAboveMin(hspectraEnergyTot,1./evts);
0812   std::cout << maxEPlot << std::endl;
0813   
0814   PlotSimple1D(canvas1DSimple, hspectraEnergyTot, -10000, maxEPlot, textSizeRel, Form("%s/EnergyTot.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, Form("#LT E_{Tot} #GT = %.1f #pm %.1f (mip/tile eq.)", hspectraEnergyTot->GetMean(), hspectraEnergyTot->GetRMS() ));
0815 
0816   
0817   hspectraEnergyTotMuon->Rebin(4);
0818   hspectraEnergyTotMuon->Scale(1./evts);
0819   hspectraEnergyTotNonMuon->Rebin(4);
0820   hspectraEnergyTotNonMuon->Scale(1./evts);
0821 
0822   TF1* fitGaussETotNonMuon = new TF1("fitEtot","gaus",hspectraEnergyTotNonMuon->GetMean()-1.5*hspectraEnergyTotNonMuon->GetRMS(),hspectraEnergyTotNonMuon->GetMean()+1.5*hspectraEnergyTotNonMuon->GetRMS());
0823   hspectraEnergyTotNonMuon->Fit(fitGaussETotNonMuon,"QRMEN0S");
0824   std::cout<< "Non muon trigg: "<<fitGaussETotNonMuon->GetParameter(1) << "\t" << fitGaussETotNonMuon->GetParameter(2) << std::endl;
0825   hspectraEnergyTotNonMuon->Fit(fitGaussETotNonMuon,"QRMEN0S","",fitGaussETotNonMuon->GetParameter(1)-1.*fitGaussETotNonMuon->GetParameter(2), fitGaussETotNonMuon->GetParameter(1)+1.*fitGaussETotNonMuon->GetParameter(2));
0826   std::cout<< "Non muon trigg refit: "<<fitGaussETotNonMuon->GetParameter(1) << "\t" << fitGaussETotNonMuon->GetParameter(2) << std::endl;
0827   
0828   
0829   PlotContamination1D(canvas1DSimple, hspectraEnergyTot,hspectraEnergyTotMuon, hspectraEnergyTotNonMuon, -10000, maxEPlot, textSizeRel, Form("%s/EnergyTotSplit.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, Form("#LT E_{Tot,non muon} #GT = %.1f#pm %.1f (mip/tile eq.)", hspectraEnergyTotNonMuon->GetMean(), hspectraEnergyTotNonMuon->GetRMS() ));
0830 
0831   PlotSimpleWithFit1D(canvas1DSimple, hspectraEnergyTotNonMuon, fitGaussETotNonMuon , -10000, maxEPlot, textSizeRel, Form("%s/EnergyTotNonMuonWithFit.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, Form("#LT E_{Tot,non muon} #GT = %.1f #pm %.1f (mip/tile eq.)", fitGaussETotNonMuon->GetParameter(1), fitGaussETotNonMuon->GetParameter(2) ));
0832 
0833   hspectraNCells->Scale(1./evts);
0834   hspectraNCells->GetYaxis()->SetTitle("counts/event");
0835   PlotSimple1D(canvas1DSimple, hspectraNCells, -10000, -10000, textSizeRel, Form("%s/NCells.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, Form("#LT N_{Cells} #GT = %.1f",hspectraNCells->GetMean() ));
0836   hspectraNCellsMuon->Scale(1./evts);
0837   hspectraNCellsNonMuon->Scale(1./evts);  
0838   PlotContamination1D( canvas1DSimple, hspectraNCells, hspectraNCellsMuon, hspectraNCellsNonMuon, -10000, -10000, textSizeRel, Form("%s/NCellsSplit.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, Form("#LT N_{Cells,non muon} #GT = %.1f",hspectraNCellsNonMuon->GetMean() ));
0839   
0840   // canvas1DSimple->SetLogy();
0841   PlotLayerOverlay(canvas1DSimple, histELayer_All, evts*10, maxXLAll ,GetAverageLayer(layersMeanAll), GetMaxLayer(layersMeanAll), textSizeRel, Form("%s/ELayerOverlay_AllTrigg.%s", outputDirPlots.Data(), plotSuffix.Data()),it->second, 1, "All triggers");
0842   PlotLayerOverlay(canvas1DSimple, histELayer_Muon, evtsMuon*10, maxXLMuon, GetAverageLayer(layersMeanMuon), GetMaxLayer(layersMeanMuon),textSizeRel, Form("%s/ELayerOverlay_MuonTrigg.%s", outputDirPlots.Data(), plotSuffix.Data()),it->second, 1, "Muon triggers");
0843   PlotLayerOverlay(canvas1DSimple, histELayer_WoMuon, (evts-evtsMuon)*10, maxXLAll, GetAverageLayer(layersMeanWOMuon), GetMaxLayer(layersMeanWOMuon),textSizeRel, Form("%s/ELayerOverlay_NonMuonTrigg.%s", outputDirPlots.Data(), plotSuffix.Data()),it->second, 1, "Non muon triggers");
0844   
0845   if (species == 0 || species == 2 ){
0846     PlotLayerOverlay(canvas1DSimple, histELayer_Muon, evtsMuon*10, maxXLMuon, GetAverageLayer(layersMeanMuon), GetMaxLayer(layersMeanMuon),textSizeRel, Form("%s/ELayerOverlay_Propaganda.%s", outputDirPlots.Data(), plotSuffix.Data()),it->second, 1, "", false);
0847   } else {
0848     PlotLayerOverlay(canvas1DSimple, histELayer_WoMuon, (evts-evtsMuon)*10, maxXLAll, GetAverageLayer(layersMeanWOMuon), GetMaxLayer(layersMeanWOMuon),textSizeRel, Form("%s/ELayerOverlay_Propaganda.%s", outputDirPlots.Data(), plotSuffix.Data()),it->second, 1, "", false);
0849   }
0850 
0851   
0852   PlotLayerOverlay(canvas1DSimple, histXLayer_All, evts*100, maxXPos-2. ,GetAverageLayer(layersMeanAll), GetMaxLayer(layersMeanAll), textSizeRel, Form("%s/XPosLayerOverlay_AllTrigg.%s", outputDirPlots.Data(), plotSuffix.Data()),it->second, 1, "All triggers");
0853   PlotLayerOverlay(canvas1DSimple, histXLayer_Muon, evtsMuon*100, maxXPos-2., GetAverageLayer(layersMeanMuon), GetMaxLayer(layersMeanMuon),textSizeRel, Form("%s/XPosLayerOverlay_MuonTrigg.%s", outputDirPlots.Data(), plotSuffix.Data()),it->second, 1, "Muon triggers");
0854   PlotLayerOverlay(canvas1DSimple, histXLayer_WoMuon, (evts-evtsMuon)*100, maxXPos-2., GetAverageLayer(layersMeanWOMuon), GetMaxLayer(layersMeanWOMuon),textSizeRel, Form("%s/XPosLayerOverlay_NonMuonTrigg.%s", outputDirPlots.Data(), plotSuffix.Data()),it->second, 1, "Non muon triggers");
0855 
0856   PlotLayerOverlay(canvas1DSimple, histYLayer_All, evts*100, maxYPos-2. ,GetAverageLayer(layersMeanAll), GetMaxLayer(layersMeanAll), textSizeRel, Form("%s/YPosLayerOverlay_AllTrigg.%s", outputDirPlots.Data(), plotSuffix.Data()),it->second, 1, "All triggers");
0857   PlotLayerOverlay(canvas1DSimple, histYLayer_Muon, evtsMuon*100, maxYPos-2., GetAverageLayer(layersMeanMuon), GetMaxLayer(layersMeanMuon),textSizeRel, Form("%s/YPosLayerOverlay_MuonTrigg.%s", outputDirPlots.Data(), plotSuffix.Data()),it->second, 1, "Muon triggers");
0858   PlotLayerOverlay(canvas1DSimple, histYLayer_WoMuon, (evts-evtsMuon)*100, maxYPos-2., GetAverageLayer(layersMeanWOMuon), GetMaxLayer(layersMeanWOMuon),textSizeRel, Form("%s/YPosLayerOverlay_NonMuonTrigg.%s", outputDirPlots.Data(), plotSuffix.Data()),it->second, 1, "Non muon triggers");
0859 
0860   PlotLayerOverlay(canvas1DSimple, histNCellsLayer_All, evts*100, maxCellsPerLayer-0.5 ,GetAverageLayer(layersMeanAll), GetMaxLayer(layersMeanAll), textSizeRel, Form("%s/NCellsLayerOverlay_AllTrigg.%s", outputDirPlots.Data(), plotSuffix.Data()),it->second, 1, "All triggers");
0861   PlotLayerOverlay(canvas1DSimple, histNCellsLayer_Muon, evtsMuon*100, maxCellsPerLayer-0.5, GetAverageLayer(layersMeanMuon), GetMaxLayer(layersMeanMuon),textSizeRel, Form("%s/NCellsLayerOverlay_MuonTrigg.%s", outputDirPlots.Data(), plotSuffix.Data()),it->second, 1, "Muon triggers");
0862   PlotLayerOverlay(canvas1DSimple, histNCellsLayer_WoMuon, (evts-evtsMuon)*100, maxCellsPerLayer-0.5, GetAverageLayer(layersMeanWOMuon), GetMaxLayer(layersMeanWOMuon),textSizeRel, Form("%s/NCellsLayerOverlay_NonMuonTrigg.%s", outputDirPlots.Data(), plotSuffix.Data()),it->second, 1, "Non muon triggers");
0863 
0864   if (ExtPlot > 0){
0865     gSystem->Exec("mkdir -p "+outputDirPlots+"/detailed");    
0866     calib.PrintGlobalInfo();  
0867     Double_t maxHG = ReturnMipPlotRangeDepVov(calib.GetVov(),true, ReadOut::Type::Caen);
0868     Double_t maxLG = 3800;
0869     
0870     //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++  
0871     // Single layer plotting
0872     //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++  
0873     DetConf::Type detConf = setup->GetDetectorConfig();
0874     MultiCanvas panelPlot(detConf, "QAData");
0875     bool init1D = panelPlot.Initialize(1);
0876     MultiCanvas panelPlot2D(detConf, "QAData2D");
0877     bool init2D = panelPlot2D.Initialize(2);
0878 
0879     std::cout << "plotting single layers" << std::endl;
0880     panelPlot.PlotMipWithFits( hSpectra, hSpectraTrigg, 1, -100, 3800, 1.2, 
0881                                 Form("%s/detailed/LocTriggerMip_HG_%.0f_%.0f" ,outputDirPlots.Data(),timemin, timemax), plotSuffix.Data(), it->second, &calib);
0882     panelPlot.PlotMipWithFits( hSpectra, hSpectraTrigg, 1, -100, maxHG, 1.2, 
0883                                 Form("%s/detailed/LocTriggerMip_Zoomed_HG_%.0f_%.0f" ,outputDirPlots.Data(),timemin, timemax), plotSuffix.Data(), it->second, &calib);    
0884     panelPlot.PlotTriggerPrim( hSpectra, averageScale, 0.8, 2., 0, 6000, 1.2,
0885                               Form("%s/detailed/All_TriggPrimitive_%.0f_%.0f" ,outputDirPlots.Data(),timemin, timemax), plotSuffix.Data(), it->second, &calib);
0886     if (ExtPlot > 1){
0887       panelPlot.PlotMipWithFits( hSpectra, hSpectraTrigg, 0, -20, maxLG, 1.2, 
0888                                 Form("%s/detailed/LocTriggerMip_LG_%.0f_%.0f" ,outputDirPlots.Data(),timemin, timemax),plotSuffix.Data(), it->second, &calib);
0889       panelPlot.PlotTriggerPrim( hSpectraTrigg, averageScale, 0.8, 2., 0, maxHG*2, 1.2,
0890                               Form("%s/detailed/LocalMuon_TriggPrimitive_%.0f_%.0f" ,outputDirPlots.Data(),timemin, timemax), plotSuffix.Data(), it->second, &calib);
0891       panelPlot2D.PlotCorrWithFits( hSpectra, 0, -20, 800, 0., 3400,
0892                                     Form("%s/detailed/LGHG_Corr_%.0f_%.0f" ,outputDirPlots.Data(),timemin, timemax), plotSuffix.Data(), it->second, &calib);
0893     }
0894     std::cout << "done plotting" << std::endl;
0895   }
0896   RootOutputHist->Close();
0897   return true;
0898 } // end QAData()
0899 
0900 
0901 //***********************************************************************************************
0902 //*********************** Simple QA routine ***************************************************
0903 //***********************************************************************************************
0904 bool DataAnalysis::SimpleQAData(void){
0905   std::cout<<"QA data"<<std::endl;
0906 
0907   std::map<int,RunInfo> ri=readRunInfosFromFile(RunListInputName.Data(),debug,0);
0908   TdataIn->GetEntry(0);
0909   Int_t runNr = event.GetRunNumber();
0910   // Get run info object
0911   std::map<int,RunInfo>::iterator it=ri.find(runNr);
0912 
0913   int species = -1;
0914   species = GetSpeciesIntFromRunInfo(it->second);
0915   float beamE = it->second.energy;
0916   std::cout << "Beam energy:" << beamE << std::endl;
0917   if (species == -1){
0918       std::cout << "WARNING: species unknown: " << it->second.species.Data() << "  aborting"<< std::endl;
0919       return false;
0920   }
0921   
0922   // create HG and LG histo's per channel
0923   TH1D* hDeltaTime = new TH1D("hDeltaTime", "Time Difference between Events; Delta Time (#mus); Counts", 2000, 0, 100000);
0924   TH2D* hspectraHGCorrvsCellID      = new TH2D( "hspectraHGCorr_vsCellID","ADC spectrum High Gain corrected vs CellID; cell ID; ADC_{HG} (arb. units)  ; counts ",
0925                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5, 4000,-200,3800);
0926   hspectraHGCorrvsCellID->SetDirectory(0);
0927   TH2D* hspectraLGCorrvsCellID      = new TH2D( "hspectraLGCorr_vsCellID","ADC spectrum Low Gain corrected vs CellID; cell ID; ADC_{LG} (arb. units)  ; counts  ",
0928                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5, 4000,-200,3800);
0929   hspectraLGCorrvsCellID->SetDirectory(0);
0930   TH2D* hspectraEnergyvsCellID  = new TH2D( "hspectraEnergy_vsCellID","Energy vs CellID; cell ID; E (mip eq./tile)    ; counts",
0931                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5, 8000,0,250);
0932   hspectraEnergyvsCellID->SetDirectory(0);
0933   
0934   std::map<int,TileSpectra> hSpectra;
0935   std::map<int, TileSpectra>::iterator ithSpectra;
0936 
0937   RootOutputHist->mkdir("IndividualCells");
0938   
0939   std::cout << "starting to run QA"<< std::endl;
0940   TcalibIn->GetEntry(0);
0941   int actCh1st = 0;
0942   double averageScale = calib.GetAverageScaleHigh(actCh1st);
0943   std::cout << "average HG mip: " << averageScale << "\t active ch: "<< actCh1st<< std::endl;
0944   
0945   int evts=TdataIn->GetEntries();
0946   double last_time = 0;
0947   double DeltaTime = 1000000000;
0948   
0949   
0950   if(percentmax!=100 || percentmin!=0){
0951     std::vector<double> deltat;
0952     deltat.reserve(evts);
0953     // Percentile calculations
0954     for (int i=0; i<evts; i++){
0955       TdataIn->GetEntry(i);
0956       
0957       double current_time = event.GetTimeStamp();
0958       if(last_time != 0){
0959         DeltaTime = current_time - last_time;
0960       }
0961       if (DeltaTime != 1000000000){
0962         deltat.push_back(DeltaTime);
0963       }
0964       last_time = current_time;
0965     }
0966     std::sort(deltat.begin(), deltat.end());
0967     int eventmax = static_cast<int>(std::round(percentmax*evts/100.0));
0968     int eventmin = static_cast<int>(std::round(percentmin*evts/100.0));
0969     std::cout << "eventmax: " << eventmax <<std::endl;
0970     std::cout <<"eventmin: " << eventmin << std::endl;
0971     if (percentmax != 100) timemax = deltat[eventmax-2];
0972     if (percentmin != 0) timemin = deltat[eventmin-1];
0973     std::cout << "timemax: " << timemax <<std::endl;
0974     std::cout <<"timemin: " << timemin << std::endl;
0975   }
0976 
0977   int rejected = 0;
0978   std::cerr<<"Run Number: " << runNr<< std::endl;
0979   std::cerr<< "Number of total evts: " << evts <<std::endl;
0980   // Event loop
0981   for(int i=0; i<evts; i++){
0982     TdataIn->GetEntry(i);
0983     
0984     // Find and fill delta time
0985     double current_time = event.GetTimeStamp();
0986     if(last_time != 0){
0987       DeltaTime = current_time - last_time;
0988       if (debug == 1000){
0989         std::cerr<< "current timestamp: "<< current_time<<std::endl;
0990       }
0991     }
0992     if (DeltaTime == 0){
0993       if(debug == 1001){
0994         std::cerr<< "Run Number: " << runNr <<"      Event Number: "<< i << std::endl;
0995         std::cerr<< "Previous Timestamp (us): "<< last_time <<"      Current Timestamp: "<< current_time<<std::endl;
0996       }
0997     }
0998         //cut based on DeltaTime range
0999     if(DeltaTime < timemin || DeltaTime > timemax){
1000       rejected++;
1001       last_time = current_time;
1002       continue;
1003     }
1004     hDeltaTime->Fill(DeltaTime);
1005     last_time = current_time;
1006     
1007     if (i%5000 == 0 && debug > 0) std::cout << "Reading " <<  i << " / " << evts << " events" << std::endl;
1008     for(int j=0; j<event.GetNTiles(); j++){
1009       Caen* aTile=(Caen*)event.GetTile(j);
1010       double corrHG = aTile->GetADCHigh()-calib.GetPedestalMeanH(aTile->GetCellID());
1011       double corrLG = aTile->GetADCLow()-calib.GetPedestalMeanL(aTile->GetCellID());
1012       hspectraHGCorrvsCellID->Fill(aTile->GetCellID(), corrHG);
1013       hspectraLGCorrvsCellID->Fill(aTile->GetCellID(), corrLG);
1014       if(aTile->GetE()!=0 ){ 
1015         hspectraEnergyvsCellID->Fill(aTile->GetCellID(), aTile->GetE());
1016       }      
1017       ithSpectra=hSpectra.find(aTile->GetCellID());
1018       if (ithSpectra!=hSpectra.end()){
1019         ithSpectra->second.FillSpectraCAEN(corrLG,corrHG);
1020         ithSpectra->second.FillTrigger(aTile->GetLocalTriggerPrimitive());;
1021         ithSpectra->second.FillCorrCAEN(corrLG,corrHG);
1022       } else {
1023         RootOutputHist->cd("IndividualCells");
1024         hSpectra[aTile->GetCellID()]=TileSpectra("AllTriggers",aTile->GetCellID(),calib.GetTileCalib(aTile->GetCellID()),event.GetROtype(),debug);
1025         hSpectra[aTile->GetCellID()].FillSpectraCAEN(corrLG,corrHG);;
1026         hSpectra[aTile->GetCellID()].FillTrigger(aTile->GetLocalTriggerPrimitive());;
1027         hSpectra[aTile->GetCellID()].FillCorrCAEN(corrLG,corrHG);
1028       }
1029     }
1030   }
1031   if (debug == 1000){
1032     std::cerr<< "Number rejected: " << rejected <<std::endl;
1033     std::cerr<< "Number accepted: " << evts - rejected << std::endl;
1034   }
1035   if (IsCalibSaveToFile()){
1036     TString fileCalibPrint = RootOutputName;
1037     fileCalibPrint         = fileCalibPrint.ReplaceAll(".root","_calib.txt");
1038     calib.PrintCalibToFile(fileCalibPrint);
1039   }
1040   
1041   RootInput->Close();      
1042   
1043   RootOutputHist->cd();
1044 
1045   hspectraHGCorrvsCellID->Write();
1046   hspectraLGCorrvsCellID->Write();
1047   hspectraEnergyvsCellID->Write();
1048 
1049   hDeltaTime->Write();
1050 
1051   RootOutputHist->cd("IndividualCells");
1052   for(ithSpectra=hSpectra.begin(); ithSpectra!=hSpectra.end(); ++ithSpectra){
1053     ithSpectra->second.Write(false);
1054   }
1055   
1056   //**********************************************************************
1057   //********************* Plotting ***************************************
1058   //**********************************************************************  
1059   TString outputDirPlots = GetPlotOutputDir();
1060   gSystem->Exec("mkdir -p "+outputDirPlots);
1061   
1062   //**********************************************************************
1063   // Create canvases for channel overview plotting
1064   //**********************************************************************
1065   Double_t textSizeRel = 0.035;
1066   StyleSettingsBasics("pdf");
1067   SetPlotStyle();
1068   
1069   TCanvas* canvasDeltaTime = new TCanvas("canvasDeltaTime","",0,0,1450,1300);  // gives the page size
1070   DefaultCanvasSettings( canvasDeltaTime, 0.08, 0.07, 0.01, 0.07);
1071   canvasDeltaTime->SetLogy(1);
1072   
1073   if(DeltaTimePlot>0){
1074     PlotSimple1D( canvasDeltaTime, hDeltaTime, -10000, timemax, textSizeRel, Form("%s/deltaTime_%.0f_%.0f.%s", outputDirPlots.Data(),timemin, timemax, plotSuffix.Data()), it->second, 1);
1075   } 
1076   
1077   TCanvas* canvas2DCorr = new TCanvas("canvasCorrPlots","",0,0,1450,1300);  // gives the page size
1078   DefaultCanvasSettings( canvas2DCorr, 0.08, 0.13, 0.045, 0.07);
1079   canvas2DCorr->SetLogz(1);
1080   TCanvas* canvas2DCorrWOLine = new TCanvas("canvasCorrPlotsWoLine","",0,0,1450,1300);  // gives the page size
1081   DefaultCanvasSettings( canvas2DCorrWOLine, 0.08, 0.13, 0.01, 0.07);
1082   canvas2DCorrWOLine->SetLogz(1);
1083   PlotSimple2D( canvas2DCorr, hspectraHGCorrvsCellID, -10000, -10000, textSizeRel, Form("%s/HGCorr_%.0f_%.0f.%s", outputDirPlots.Data(),timemin, timemax, plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
1084   PlotSimple2D( canvas2DCorr, hspectraLGCorrvsCellID, -10000, -10000, textSizeRel, Form("%s/LGCorr_%.0f_%.0f.%s", outputDirPlots.Data(),timemin, timemax, plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
1085   PlotSimple2D( canvas2DCorr, hspectraEnergyvsCellID, -10000, -10000, textSizeRel, Form("%s/EnergyVsCellID_%.0f_%.0f.%s", outputDirPlots.Data(),timemin, timemax, plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
1086     
1087   if (ExtPlot > 0){
1088     gSystem->Exec("mkdir -p "+outputDirPlots+"/detailed");
1089     //***********************************************************************************************************
1090     //********************************* 8 Panel overview plot  **************************************************
1091     //***********************************************************************************************************
1092     //*****************************************************************
1093       // Test beam geometry (beam coming from viewer)
1094       //===========================================================
1095       //||    8 (4)    ||    7 (5)   ||    6 (6)   ||    5 (7)   ||  row 0
1096       //===========================================================
1097       //||    1 (0)    ||    2 (1)   ||    3 (2)   ||    4 (3)   ||  row 1
1098       //===========================================================
1099       //    col 0     col 1       col 2     col  3
1100       // rebuild pad geom in similar way (numbering -1)
1101     //*****************************************************************
1102     TCanvas* canvas8Panel;
1103     TPad* pad8Panel[8];
1104     Double_t topRCornerX[8];
1105     Double_t topRCornerY[8];
1106     Int_t textSizePixel = 30;
1107     Double_t relSize8P[8];
1108     CreateCanvasAndPadsFor8PannelTBPlot(canvas8Panel, pad8Panel,  topRCornerX, topRCornerY, relSize8P, textSizePixel);
1109   
1110     TCanvas* canvas8PanelProf;
1111     TPad* pad8PanelProf[8];
1112     Double_t topRCornerXProf[8];
1113     Double_t topRCornerYProf[8];
1114     Double_t relSize8PProf[8];
1115     CreateCanvasAndPadsFor8PannelTBPlot(canvas8PanelProf, pad8PanelProf,  topRCornerXProf, topRCornerYProf, relSize8PProf, textSizePixel, 0.045, "Prof", false);
1116 
1117     
1118     calib.PrintGlobalInfo();  
1119     Double_t maxHG = ReturnMipPlotRangeDepVov(calib.GetVov(),true, ReadOut::Type::Caen);
1120     Double_t maxLG = 3800;
1121     
1122         //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++  
1123     // Single layer plotting
1124     //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++  
1125     DetConf::Type detConf = setup->GetDetectorConfig();
1126     MultiCanvas panelPlot(detConf, "QAData");
1127     bool init1D = panelPlot.Initialize(1);
1128     MultiCanvas panelPlot2D(detConf, "QAData2D");
1129     bool init2D = panelPlot2D.Initialize(2);
1130 
1131     std::cout << "plotting single layers" << std::endl;
1132     panelPlot.PlotNoiseWithFits( hSpectra,  0, 0, maxHG, 1.2, 
1133                                 Form("%s/detailed/Spectra_HG_%.0f_%.0f" ,outputDirPlots.Data(),timemin, timemax), plotSuffix.Data(), it->second, &calib);
1134     panelPlot.PlotNoiseWithFits( hSpectra,  1, 0, maxLG, 1.2, 
1135                                 Form("%s/detailed/Spectra_LG_%.0f_%.0f" ,outputDirPlots.Data(),timemin, timemax), plotSuffix.Data(), it->second, &calib);
1136     panelPlot2D.PlotCorrWithFits( hSpectra, 0, -20, 340, 0., 3800,
1137                                   Form("%s/detailed/LGHG_Corr_%.0f_%.0f" ,outputDirPlots.Data(),timemin, timemax), plotSuffix.Data(), it->second, &calib);
1138     std::cout << "done plotting" << std::endl;
1139   }
1140   RootOutputHist->Close();
1141   return true;
1142 }
1143 
1144 //***********************************************************************************************
1145 //*********************** Create output files ***************************************************
1146 //***********************************************************************************************
1147 bool DataAnalysis::CreateOutputRootFile(void){
1148   if(Overwrite){
1149     RootOutput=new TFile(RootOutputName.Data(),"RECREATE");
1150   } else{
1151     RootOutput = new TFile(RootOutputName.Data(),"CREATE");
1152   }
1153   if(RootOutput->IsZombie()){
1154     std::cout<<"Error opening '"<<RootOutput<<"'no reachable path? Exist without force mode to overwrite?..."<<std::endl;
1155     return false;
1156   }
1157   return true;
1158 }
1159 
1160 //***********************************************************************************************
1161 //*********************** Create output files ***************************************************
1162 //***********************************************************************************************
1163 bool DataAnalysis::CreateOutputRootFileHist(void){
1164   if(Overwrite){
1165     RootOutputHist=new TFile(RootOutputNameHist.Data(),"RECREATE");
1166   } else{
1167     RootOutputHist = new TFile(RootOutputNameHist.Data(),"CREATE");
1168   }
1169   if(RootOutputHist->IsZombie()){
1170     std::cout<<"Error opening '"<<RootOutputHist<<"'no reachable path? Exist without force mode to overwrite?..."<<std::endl;
1171     return false;
1172   }
1173   return true;
1174 }
[ 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!