Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-12-16 09:28:18 

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