File indexing completed on 2025-07-13 07:54:15
0001 #include "DataAnalysis.h"
0002 #include <vector>
0003 #include "TROOT.h"
0004
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
0022
0023 bool DataAnalysis::CheckAndOpenIO(void){
0024 int matchingbranch;
0025
0026
0027 std::cout<<"Input name set to: '"<<RootInputName.Data() <<std::endl;
0028 if(!RootInputName.IsNull()){
0029
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
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
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
0056
0057
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
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
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
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
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
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
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
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
0174 TH2D* hspectraEnergyTotvsNCells = nullptr;
0175 TH2D* hspectraEnergyTotvsNCellsMuon = nullptr;
0176 TH2D* hspectraEnergyVsLayer = nullptr;
0177 TH2D* hspectraEnergyVsLayerMuon = nullptr;
0178
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
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
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
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
0342
0343
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
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
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
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
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
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
0584
0585 TString outputDirPlots = GetPlotOutputDir();
0586 gSystem->Exec("mkdir -p "+outputDirPlots);
0587
0588
0589
0590
0591 Double_t textSizeRel = 0.035;
0592 StyleSettingsBasics("pdf");
0593 SetPlotStyle();
0594
0595 TCanvas* canvas2DCorr = new TCanvas("canvasCorrPlots","",0,0,1450,1300);
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);
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);
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
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
0691
0692
0693
0694
0695
0696
0697
0698
0699
0700
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
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
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
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
0846
0847 TString outputDirPlots = GetPlotOutputDir();
0848 gSystem->Exec("mkdir -p "+outputDirPlots);
0849
0850
0851
0852
0853 Double_t textSizeRel = 0.035;
0854 StyleSettingsBasics("pdf");
0855 SetPlotStyle();
0856
0857 TCanvas* canvas2DCorr = new TCanvas("canvasCorrPlots","",0,0,1450,1300);
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);
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
0871
0872
0873
0874
0875
0876
0877
0878
0879
0880
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
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
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 }