Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-07-13 07:54:15

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