Back to home page

EIC code displayed by LXR
 
 

    

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

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   // copy full calibration to different file and calculate energy
0124   if(RunQA){
0125     status=QAData();
0126   }  
0127   return status;
0128 }
0129 
0130 
0131 
0132 //***********************************************************************************************
0133 //*********************** Calibration routine ***************************************************
0134 //***********************************************************************************************
0135 bool DataAnalysis::QAData(void){
0136   std::cout<<"QA data"<<std::endl;
0137 
0138   std::map<int,RunInfo> ri=readRunInfosFromFile(RunListInputName.Data(),debug,0);
0139   
0140   // create HG and LG histo's per channel
0141   TH2D* hspectraHGCorrvsCellID      = new TH2D( "hspectraHGCorr_vsCellID","ADC spectrum High Gain corrected vs CellID; cell ID; ADC_{HG} (arb. units)  ; counts ",
0142                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5, 4000,-200,3800);
0143   hspectraHGCorrvsCellID->SetDirectory(0);
0144   TH2D* hspectraLGCorrvsCellID      = new TH2D( "hspectraLGCorr_vsCellID","ADC spectrum Low Gain corrected vs CellID; cell ID; ADC_{LG} (arb. units)  ; counts  ",
0145                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5, 4000,-200,3800);
0146   hspectraLGCorrvsCellID->SetDirectory(0);
0147   TH2D* hspectraHGvsCellID      = new TH2D( "hspectraHG_vsCellID","ADC spectrum High Gain vs CellID; cell ID; ADC_{HG} (arb. units)   ; counts ",
0148                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5, 4000,0,4000);
0149   hspectraHGvsCellID->SetDirectory(0);
0150   TH2D* hspectraLGvsCellID      = new TH2D( "hspectraLG_vsCellID","ADC spectrum Low Gain vs CellID; cell ID; ADC_{LG} (arb. units)  ; counts",
0151                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5, 4000,0,4000);
0152   hspectraLGvsCellID->SetDirectory(0);
0153   TH2D* hspectraEnergyvsCellID  = new TH2D( "hspectraEnergy_vsCellID","Energy vs CellID; cell ID; E (mip eq./tile)    ; counts",
0154                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5, 6000,0,200);
0155   hspectraEnergyvsCellID->SetDirectory(0);
0156   TH2D* hspectraEnergyTotvsNCells  = new TH2D( "hspectraTotEnergy_vsNCells","Energy vs CellID; N_{Cells}; E_{tot} (mip eq./tile) ; counts",
0157                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5, 6000,0,1000);
0158   hspectraEnergyTotvsNCells->SetDirectory(0);
0159   TH2D* hspectraEnergyTotvsNCellsMuon  = new TH2D( "hspectraTotEnergy_vsNCellsMuon","Energy vs CellID; N_{Cells}; E_{tot} (mip eq./tile) ; counts",
0160                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5, 6000,0,1000);
0161   hspectraEnergyTotvsNCellsMuon->SetDirectory(0);
0162 
0163 
0164   TH2D* hspectraEnergyVsLayer  = new TH2D( "hspectraTotLayerEnergy_vsLayer","Energy in layer vs Layer; Layer; E_{layer} (mip eq./tile) ; counts",
0165                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, 6000,0,1000);
0166   hspectraEnergyVsLayer->SetDirectory(0);
0167   TH2D* hspectraEnergyVsLayerMuon  = new TH2D( "hspectraTotLayerEnergy_vsLayerMuon","Energy in layer vs Layer Muon triggers; Layer; E_{layer} (mip eq./tile) ; counts",
0168                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, 6000,0,1000);
0169   hspectraEnergyVsLayerMuon->SetDirectory(0);
0170   TH2D* hAverageXVsLayer  = new TH2D( "hAverageX_vsLayer","Av. X pos in layer vs Layer; Layer; X_{pos} (cm) ; counts",
0171                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, 100,-10,10);
0172   hAverageXVsLayer->SetDirectory(0);
0173   TH2D* hAverageYVsLayer  = new TH2D( "hAverageX_vsLayer","Av. Y pos in layer vs Layer; Layer; Y_{pos} (cm) ; counts",
0174                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, 50,-5,5);
0175   hAverageYVsLayer->SetDirectory(0);
0176   TH2D* hNCellsVsLayer  = new TH2D( "hnCells_vsLayer","NCells in layer vs Layer; Layer; N_{cells,layer} ; counts",
0177                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, 9,-0.5,8.5);
0178   
0179   
0180   
0181   std::map<int,TileSpectra> hSpectra;
0182   std::map<int,TileSpectra> hSpectraTrigg;
0183   std::map<int, TileSpectra>::iterator ithSpectra;
0184   std::map<int, TileSpectra>::iterator ithSpectraTrigg;
0185 
0186   RootOutputHist->mkdir("IndividualCells");
0187   RootOutputHist->mkdir("IndividualCellsTrigg");
0188   
0189   Int_t runNr = -1;
0190   std::cout << "starting to run QA"<< std::endl;
0191   TcalibIn->GetEntry(0);
0192   int actCh1st = 0;
0193   double averageScale = calib.GetAverageScaleHigh(actCh1st);
0194   std::cout << "average HG mip: " << averageScale << "\t active ch: "<< actCh1st<< std::endl;
0195   
0196   int evts=TdataIn->GetEntries();
0197   int evtsMuon= 0;
0198   for(int i=0; i<evts; i++){
0199     TdataIn->GetEntry(i);
0200     if (i%5000 == 0 && debug > 0) std::cout << "Reading " <<  i << " / " << evts << " events" << std::endl;
0201     if (i == 0)runNr = event.GetRunNumber();
0202     double Etot = 0;
0203     int nCells  = 0;
0204     bool muontrigg = false;
0205     int locMuon = 0;
0206     std::map<int,Layer> layers;
0207     std::map<int, Layer>::iterator ithLayer;
0208     
0209     for(int j=0; j<event.GetNTiles(); j++){
0210       Caen* aTile=(Caen*)event.GetTile(j);
0211       double corrHG = aTile->GetADCHigh()-calib.GetPedestalMeanH(aTile->GetCellID());
0212       double corrLG = aTile->GetADCLow()-calib.GetPedestalMeanL(aTile->GetCellID());
0213       hspectraHGvsCellID->Fill(aTile->GetCellID(), aTile->GetADCHigh());
0214       hspectraLGvsCellID->Fill(aTile->GetCellID(), aTile->GetADCLow());
0215       hspectraHGCorrvsCellID->Fill(aTile->GetCellID(), corrHG);
0216       hspectraLGCorrvsCellID->Fill(aTile->GetCellID(), corrLG);
0217       if(aTile->GetE()!=0){ 
0218         nCells++;
0219         int currLayer = setup->GetLayer(aTile->GetCellID());
0220         ithLayer=layers.find(currLayer);
0221         if(ithLayer!=layers.end()){
0222           ithLayer->second.nCells+=1;
0223           ithLayer->second.energy+=aTile->GetE();
0224           ithLayer->second.avX+=setup->GetX(aTile->GetCellID())*aTile->GetE();
0225           ithLayer->second.avY+=setup->GetY(aTile->GetCellID())*aTile->GetE();
0226         } else {
0227           layers[currLayer]=Layer();
0228           layers[currLayer].nCells+=1;
0229           layers[currLayer].energy+=aTile->GetE();
0230           layers[currLayer].avX+=setup->GetX(aTile->GetCellID())*aTile->GetE();
0231           layers[currLayer].avY+=setup->GetY(aTile->GetCellID())*aTile->GetE();
0232         }
0233       }
0234       // read current tile
0235       ithSpectraTrigg=hSpectraTrigg.find(aTile->GetCellID());
0236       if (aTile->GetLocalTriggerBit() == 1){
0237         if(ithSpectraTrigg!=hSpectraTrigg.end()){
0238           ithSpectraTrigg->second.FillTrigger(aTile->GetLocalTriggerPrimitive());
0239           ithSpectraTrigg->second.FillSpectra(corrLG,corrHG);
0240         } else {
0241           RootOutputHist->cd("IndividualCellsTrigg");
0242           hSpectraTrigg[aTile->GetCellID()]=TileSpectra("mipTrigg",aTile->GetCellID(),calib.GetTileCalib(aTile->GetCellID()),debug);
0243           hSpectraTrigg[aTile->GetCellID()].FillTrigger(aTile->GetLocalTriggerPrimitive());
0244           hSpectraTrigg[aTile->GetCellID()].FillSpectra(corrLG,corrHG);
0245         }
0246         locMuon++;        
0247       }      
0248       
0249       ithSpectra=hSpectra.find(aTile->GetCellID());
0250       if (ithSpectra!=hSpectra.end()){
0251         ithSpectra->second.FillSpectra(corrLG,corrHG);
0252         ithSpectra->second.FillTrigger(aTile->GetLocalTriggerPrimitive());;
0253         ithSpectra->second.FillCorr(corrLG,corrHG);
0254       } else {
0255         RootOutputHist->cd("IndividualCells");
0256         hSpectra[aTile->GetCellID()]=TileSpectra("mip1st",aTile->GetCellID(),calib.GetTileCalib(aTile->GetCellID()),debug);
0257         hSpectra[aTile->GetCellID()].FillSpectra(corrLG,corrHG);;
0258         hSpectra[aTile->GetCellID()].FillTrigger(aTile->GetLocalTriggerPrimitive());;
0259         hSpectra[aTile->GetCellID()].FillCorr(corrLG,corrHG);
0260       }
0261     }
0262     if (nCells > 0) {
0263 
0264       int nLayerSingleCell = 0;
0265       for(ithLayer=layers.begin(); ithLayer!=layers.end(); ++ithLayer){
0266         if (ithLayer->second.nCells == 1) 
0267             nLayerSingleCell++;
0268       }
0269       double fracLayer1Cell = (double)nLayerSingleCell/(int)layers.size();
0270       // if ( fracLayer1Cell > 0.8){
0271       //   muontrigg = true;
0272       //   evtsMuon++;        
0273       // }
0274       double fracLocMuon = (double)locMuon/nCells;
0275       if (fracLocMuon > 0.6){ 
0276         muontrigg = true;
0277         evtsMuon++;
0278       }
0279       if (muontrigg && debug > 3){
0280           std::cout << "Muon triggered:\t" <<  fracLayer1Cell << "\t" << fracLocMuon << std::endl;
0281       }
0282       
0283       for(int l = 0; l < setup->GetNMaxLayer()+1;l++){
0284         ithLayer=layers.find(l);
0285         if (ithLayer!=layers.end()){          
0286           double avx     = ithLayer->second.avX/ithLayer->second.energy;
0287           double avy     = ithLayer->second.avY/ithLayer->second.energy;
0288           hspectraEnergyVsLayer->Fill(l,ithLayer->second.energy);
0289           if (muontrigg) hspectraEnergyVsLayerMuon->Fill(l,ithLayer->second.energy);
0290           hAverageXVsLayer->Fill(l,avx);
0291           hAverageYVsLayer->Fill(l,avy);
0292           hNCellsVsLayer->Fill(l,ithLayer->second.nCells);
0293           if (ithLayer->second.nCells == 1) 
0294             nLayerSingleCell++;
0295         } else {
0296           hspectraEnergyVsLayer->Fill(l,0.);
0297           if (muontrigg) hspectraEnergyVsLayerMuon->Fill(l,0.);
0298           hAverageXVsLayer->Fill(l,-20);
0299           hAverageYVsLayer->Fill(l,-20);
0300           hNCellsVsLayer->Fill(l,0);
0301         }
0302       }
0303 
0304       
0305       
0306       for(int j=0; j<event.GetNTiles(); j++){
0307         Caen* aTile=(Caen*)event.GetTile(j);
0308         // remove bad channels from output
0309         double energy = aTile->GetE();
0310         if(energy!=0){ 
0311           hspectraEnergyvsCellID->Fill(aTile->GetCellID(), energy);
0312           Etot=Etot+energy;
0313         } 
0314       }
0315 
0316       hspectraEnergyTotvsNCells->Fill(nCells,Etot);
0317       if(muontrigg) hspectraEnergyTotvsNCellsMuon->Fill(nCells,Etot);
0318     }
0319   }
0320   
0321   if (IsCalibSaveToFile()){
0322     TString fileCalibPrint = RootOutputName;
0323     fileCalibPrint         = fileCalibPrint.ReplaceAll(".root","_calib.txt");
0324     calib.PrintCalibToFile(fileCalibPrint);
0325   }
0326   
0327   RootInput->Close();      
0328   
0329   RootOutputHist->cd();
0330 
0331   hspectraHGvsCellID->Write();
0332   hspectraLGvsCellID->Write();
0333   hspectraHGCorrvsCellID->Write();
0334   hspectraLGCorrvsCellID->Write();
0335   hspectraEnergyvsCellID->Write();
0336   hspectraEnergyTotvsNCells->Write();
0337   hspectraEnergyTotvsNCellsMuon->Write();
0338   
0339   TH2D* hspectraEnergyTotvsNCells_WoMuon = (TH2D*)hspectraEnergyTotvsNCells->Clone("hspectraTotEnergy_vsNCells_woMuon");
0340   hspectraEnergyTotvsNCells_WoMuon->Sumw2();
0341   hspectraEnergyTotvsNCells_WoMuon->Add(hspectraEnergyTotvsNCellsMuon, -1);
0342   hspectraEnergyTotvsNCells_WoMuon->Write();
0343   
0344   hspectraEnergyTotvsNCells->Sumw2();
0345   TH1D* hspectraEnergyTot = (TH1D*)hspectraEnergyTotvsNCells->ProjectionY();
0346   hspectraEnergyTot->SetDirectory(0);
0347   TH1D* hspectraNCells = (TH1D*)hspectraEnergyTotvsNCells->ProjectionX();
0348   hspectraNCells->SetDirectory(0);
0349   hspectraEnergyTot->Write("hTotEnergy");
0350   hspectraNCells->Write("hNCells");
0351 
0352   hspectraEnergyTotvsNCellsMuon->Sumw2();
0353   TH1D* hspectraEnergyTotMuon = (TH1D*)hspectraEnergyTotvsNCellsMuon->ProjectionY();
0354   hspectraEnergyTotMuon->SetDirectory(0);
0355   TH1D* hspectraNCellsMuon = (TH1D*)hspectraEnergyTotvsNCellsMuon->ProjectionX();
0356   hspectraNCellsMuon->SetDirectory(0);
0357   hspectraEnergyTotMuon->Write("hTotEnergyMuon");
0358   hspectraNCellsMuon->Write("hNCellsMuon");
0359   
0360   TH1D* hspectraEnergyTotNonMuon = (TH1D*)hspectraEnergyTotvsNCells_WoMuon->ProjectionY();
0361   hspectraEnergyTotNonMuon->SetDirectory(0);
0362   TH1D* hspectraNCellsNonMuon = (TH1D*)hspectraEnergyTotvsNCells_WoMuon->ProjectionX();
0363   hspectraNCellsNonMuon->SetDirectory(0);
0364   hspectraEnergyTotNonMuon->Write("hTotEnergyNonMuon");
0365   hspectraNCellsNonMuon->Write("hNCellsNonMuon");
0366   
0367   hAverageXVsLayer->Write();
0368   hAverageYVsLayer->Write();
0369   hNCellsVsLayer->Write();
0370   hspectraEnergyVsLayer->Write();
0371   hspectraEnergyVsLayerMuon->Write();
0372   TH2D* hspectraEnergyVsLayer_WoMuon = (TH2D*)hspectraEnergyVsLayer->Clone("hspectraTotLayerEnergy_vsLayer_woMuon");
0373   hspectraEnergyVsLayer_WoMuon->Sumw2();
0374   hspectraEnergyVsLayer_WoMuon->Add(hspectraEnergyVsLayerMuon, -1);
0375   hspectraEnergyVsLayer_WoMuon->Write();
0376 
0377   
0378   //**********************************************************************
0379   //********************* Plotting ***************************************
0380   //**********************************************************************
0381   // Get run info object
0382   std::map<int,RunInfo>::iterator it=ri.find(runNr);
0383   
0384   TString outputDirPlots = GetPlotOutputDir();
0385   gSystem->Exec("mkdir -p "+outputDirPlots);
0386   
0387   //**********************************************************************
0388   // Create canvases for channel overview plotting
0389   //**********************************************************************
0390   Double_t textSizeRel = 0.035;
0391   StyleSettingsBasics("pdf");
0392   SetPlotStyle();
0393   
0394   TCanvas* canvas2DCorr = new TCanvas("canvasCorrPlots","",0,0,1450,1300);  // gives the page size
0395   DefaultCancasSettings( canvas2DCorr, 0.08, 0.13, 0.045, 0.07);
0396   canvas2DCorr->SetLogz(1);
0397   PlotSimple2D( canvas2DCorr, hspectraHGvsCellID, -10000, -10000, textSizeRel, Form("%s/HG.pdf", outputDirPlots.Data()), it->second, 1, kFALSE, "colz", true);
0398   PlotSimple2D( canvas2DCorr, hspectraLGvsCellID, -10000, -10000, textSizeRel, Form("%s/LG.pdf", outputDirPlots.Data()), it->second, 1, kFALSE, "colz", true);
0399   PlotSimple2D( canvas2DCorr, hspectraHGCorrvsCellID, -10000, -10000, textSizeRel, Form("%s/HGCorr.pdf", outputDirPlots.Data()), it->second, 1, kFALSE, "colz", true);
0400   PlotSimple2D( canvas2DCorr, hspectraLGCorrvsCellID, -10000, -10000, textSizeRel, Form("%s/LGCorr.pdf", outputDirPlots.Data()), it->second, 1, kFALSE, "colz", true);
0401   PlotSimple2D( canvas2DCorr, hspectraEnergyvsCellID, -10000, -10000, textSizeRel, Form("%s/EnergyVsCellID.pdf", outputDirPlots.Data()), it->second, 1, kFALSE, "colz", true);
0402   PlotSimple2D( canvas2DCorr, hspectraEnergyTotvsNCells, -10000, -10000, textSizeRel, Form("%s/EnergyTotalVsNCells.pdf", outputDirPlots.Data()), it->second, 1, kFALSE, "colz", true, Form("evts = %d",evts));
0403   PlotSimple2D( canvas2DCorr, hspectraEnergyTotvsNCellsMuon, -10000, -10000, textSizeRel, Form("%s/EnergyTotalVsNCells_MuonTrigg.pdf", outputDirPlots.Data()), it->second, 1, kFALSE, "colz", true, Form("Muon evts = %d",evtsMuon));
0404   PlotSimple2D( canvas2DCorr, hspectraEnergyTotvsNCells_WoMuon, -10000, -10000, textSizeRel, Form("%s/EnergyTotalVsNCells_WoMuonTrigg.pdf", outputDirPlots.Data()), it->second, 1, kFALSE, "colz", true, Form("Non Muon evts = %d",evts-evtsMuon));
0405 
0406 
0407   PlotSimple2D( canvas2DCorr, hspectraEnergyVsLayer, -10000, -10000, textSizeRel, Form("%s/EnergyVsLayer.pdf", outputDirPlots.Data()), it->second, 1, kFALSE, "colz", true);
0408   PlotSimple2D( canvas2DCorr, hspectraEnergyVsLayerMuon, -10000, -10000, textSizeRel, Form("%s/EnergyVsLayer_MuonTrigg.pdf", outputDirPlots.Data()), it->second, 1, kFALSE, "colz", true, "Muon triggers");
0409   PlotSimple2D( canvas2DCorr, hspectraEnergyVsLayer_WoMuon, -10000, -10000, textSizeRel, Form("%s/EnergyVsLayer_WOMuonTrigg.pdf", outputDirPlots.Data()), it->second, 1, kFALSE, "colz", true, "Non Muon triggers");
0410   PlotSimple2D( canvas2DCorr, hAverageXVsLayer, -10000, -10000, textSizeRel, Form("%s/XPosVsLayer.pdf", outputDirPlots.Data()), it->second, 1, kFALSE, "colz", true);
0411   PlotSimple2D( canvas2DCorr, hAverageYVsLayer, -10000, -10000, textSizeRel, Form("%s/YPosVsLayer.pdf", outputDirPlots.Data()), it->second, 1, kFALSE, "colz", true);
0412   PlotSimple2D( canvas2DCorr, hNCellsVsLayer, -10000, -10000, textSizeRel, Form("%s/NcellsLayerVsLayer.pdf", outputDirPlots.Data()), it->second, 1, kFALSE, "colz", true);
0413   
0414   
0415   TCanvas* canvas1DSimple = new TCanvas("canvas1DSimple","",0,0,1450,1300);  // gives the page size
0416   DefaultCancasSettings( canvas1DSimple, 0.08, 0.03, 0.03, 0.07);
0417   hspectraEnergyTot->Scale(1./evts);
0418   hspectraEnergyTot->GetYaxis()->SetTitle("counts/event");
0419   PlotSimple1D(canvas1DSimple, hspectraEnergyTot, -10000, -10000, textSizeRel, Form("%s/EnergyTot.pdf", outputDirPlots.Data()), it->second, 1, Form("#LT E_{Tot} #GT = %.1f (mip/tile eq.)",hspectraEnergyTot->GetMean() ));
0420 
0421   hspectraEnergyTotMuon->Scale(1./evts);
0422   hspectraEnergyTotNonMuon->Scale(1./evts);
0423   PlotContamination1D(canvas1DSimple, hspectraEnergyTot,hspectraEnergyTotMuon, hspectraEnergyTotNonMuon, -10000, -10000, textSizeRel, Form("%s/EnergyTotSplit.pdf", outputDirPlots.Data()), it->second, 1, Form("#LT E_{Tot,non muon} #GT = %.1f (mip/tile eq.)",hspectraEnergyTotNonMuon->GetMean() ));
0424 
0425   hspectraNCells->Scale(1./evts);
0426   hspectraNCells->GetYaxis()->SetTitle("counts/event");
0427   PlotSimple1D(canvas1DSimple, hspectraNCells, -10000, -10000, textSizeRel, Form("%s/NCells.pdf", outputDirPlots.Data()), it->second, 1, Form("#LT N_{Cells} #GT = %.1f",hspectraNCells->GetMean() ));
0428   hspectraNCellsMuon->Scale(1./evts);
0429   hspectraNCellsNonMuon->Scale(1./evts);  
0430   PlotContamination1D( canvas1DSimple, hspectraNCells, hspectraNCellsMuon, hspectraNCellsNonMuon, -10000, -10000, textSizeRel, Form("%s/NCellsSplit.pdf", outputDirPlots.Data()), it->second, 1, Form("#LT N_{Cells,non muon} #GT = %.1f",hspectraNCellsNonMuon->GetMean() ));
0431   
0432   
0433   if (ExtPlot > 0){
0434     //***********************************************************************************************************
0435     //********************************* 8 Panel overview plot  **************************************************
0436     //***********************************************************************************************************
0437     //*****************************************************************
0438       // Test beam geometry (beam coming from viewer)
0439       //===========================================================
0440       //||    8 (4)    ||    7 (5)   ||    6 (6)   ||    5 (7)   ||  row 0
0441       //===========================================================
0442       //||    1 (0)    ||    2 (1)   ||    3 (2)   ||    4 (3)   ||  row 1
0443       //===========================================================
0444       //    col 0     col 1       col 2     col  3
0445       // rebuild pad geom in similar way (numbering -1)
0446     //*****************************************************************
0447     TCanvas* canvas8Panel;
0448     TPad* pad8Panel[8];
0449     Double_t topRCornerX[8];
0450     Double_t topRCornerY[8];
0451     Int_t textSizePixel = 30;
0452     Double_t relSize8P[8];
0453     CreateCanvasAndPadsFor8PannelTBPlot(canvas8Panel, pad8Panel,  topRCornerX, topRCornerY, relSize8P, textSizePixel);
0454   
0455     TCanvas* canvas8PanelProf;
0456     TPad* pad8PanelProf[8];
0457     Double_t topRCornerXProf[8];
0458     Double_t topRCornerYProf[8];
0459     Double_t relSize8PProf[8];
0460     CreateCanvasAndPadsFor8PannelTBPlot(canvas8PanelProf, pad8PanelProf,  topRCornerXProf, topRCornerYProf, relSize8PProf, textSizePixel, 0.045, "Prof", false);
0461 
0462     
0463     calib.PrintGlobalInfo();  
0464     Double_t maxHG = ReturnMipPlotRangeDepVov(calib.GetVov(),true);
0465     Double_t maxLG = 3800;
0466     std::cout << "plotting single layers" << std::endl;
0467     for (Int_t l = 0; l < setup->GetNMaxLayer()+1; l++){    
0468       if (l%10 == 0 && l > 0 && debug > 0){
0469         std::cout << "============================== layer " <<  l << " / " << setup->GetNMaxLayer() << " layers" << std::endl;
0470       }
0471       PlotMipWithFitsFullLayer (canvas8Panel,pad8Panel, topRCornerX, topRCornerY, relSize8P, textSizePixel, 
0472                               hSpectra, hSpectraTrigg, setup, true, -100, 3800, 1.2, l, 0,
0473                               Form("%s/LocTriggerMip_HG_Layer%02d.pdf" ,outputDirPlots.Data(), l), it->second);
0474       PlotMipWithFitsFullLayer (canvas8Panel,pad8Panel, topRCornerX, topRCornerY, relSize8P, textSizePixel, 
0475                               hSpectra, hSpectraTrigg, setup, true, -100, maxHG, 1.2, l, 0,
0476                               Form("%s/LocTriggerMip_Zoomed_HG_Layer%02d.pdf" ,outputDirPlots.Data(), l), it->second);
0477       PlotTriggerPrimWithFitsFullLayer (canvas8Panel,pad8Panel, topRCornerX, topRCornerY, relSize8P, textSizePixel, 
0478                                         hSpectra, setup, averageScale, 0.8, 2.,
0479                                         0, 6000, 1.2, l, 0, Form("%s/All_TriggPrimitive_Layer%02d.pdf" ,outputDirPlots.Data(), l), it->second);
0480       if (ExtPlot > 1){
0481         PlotCorrWithFitsFullLayer(canvas8PanelProf,pad8PanelProf, topRCornerXProf, topRCornerYProf, relSize8PProf, textSizePixel, 
0482                                       hSpectra, setup, false, -20, 800, 1.2, l, 0,
0483                                       Form("%s/LGHG_Corr_Layer%02d.pdf" ,outputDirPlots.Data(), l), it->second);
0484         PlotMipWithFitsFullLayer (canvas8Panel,pad8Panel, topRCornerX, topRCornerY, relSize8P, textSizePixel, 
0485                                   hSpectra, hSpectraTrigg, setup, false, -20, maxLG, 1.2, l, 0,
0486                                   Form("%s/LocTriggerMip_LG_Layer%02d.pdf" ,outputDirPlots.Data(), l), it->second);
0487         PlotTriggerPrimWithFitsFullLayer (canvas8Panel,pad8Panel, topRCornerX, topRCornerY, relSize8P, textSizePixel, 
0488                                           hSpectraTrigg, setup, averageScale, 0.8, 2.,
0489                                           0, maxHG*2, 1.2, l, 0, Form("%s/LocalMuon_TriggPrimitive_Layer%02d.pdf" ,outputDirPlots.Data(), l), it->second);
0490       }
0491     }
0492     std::cout << "done plotting" << std::endl;
0493   }
0494   RootOutputHist->Close();
0495   return true;
0496 }
0497 
0498 //***********************************************************************************************
0499 //*********************** Create output files ***************************************************
0500 //***********************************************************************************************
0501 bool DataAnalysis::CreateOutputRootFile(void){
0502   if(Overwrite){
0503     RootOutput=new TFile(RootOutputName.Data(),"RECREATE");
0504   } else{
0505     RootOutput = new TFile(RootOutputName.Data(),"CREATE");
0506   }
0507   if(RootOutput->IsZombie()){
0508     std::cout<<"Error opening '"<<RootOutput<<"'no reachable path? Exist without force mode to overwrite?..."<<std::endl;
0509     return false;
0510   }
0511   return true;
0512 }
0513 
0514 //***********************************************************************************************
0515 //*********************** Create output files ***************************************************
0516 //***********************************************************************************************
0517 bool DataAnalysis::CreateOutputRootFileHist(void){
0518   if(Overwrite){
0519     RootOutputHist=new TFile(RootOutputNameHist.Data(),"RECREATE");
0520   } else{
0521     RootOutputHist = new TFile(RootOutputNameHist.Data(),"CREATE");
0522   }
0523   if(RootOutputHist->IsZombie()){
0524     std::cout<<"Error opening '"<<RootOutputHist<<"'no reachable path? Exist without force mode to overwrite?..."<<std::endl;
0525     return false;
0526   }
0527   return true;
0528 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!