Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-02-22 09:39:22 

0001 #include "EventDisplay.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 "TH3F.h"
0011 #include "TPaletteAxis.h"
0012 #include "TProfile.h"
0013 #include "TChain.h"
0014 #include "CommonHelperFunctions.h"
0015 #include "TileSpectra.h"
0016 #include "PlottHelper.h"
0017 #include "TRandom3.h"
0018 #include "TMath.h"
0019 #include "Math/Minimizer.h"
0020 #include "Math/MinimizerOptions.h"
0021 
0022 // ****************************************************************************
0023 // Checking and opening input and output files
0024 // ****************************************************************************
0025 bool EventDisplay::CheckAndOpenIO(void){
0026   int matchingbranch;
0027 
0028   //Need to check first input to get the setup...I do not think it is necessary
0029   std::cout<<"Input name set to: '"<<RootInputName.Data() <<std::endl;  
0030   if(!RootInputName.IsNull()){
0031     //File exist?
0032     RootInput=new TFile(RootInputName.Data(),"READ");
0033     if(RootInput->IsZombie()){
0034       std::cout<<"Error opening '"<<RootInputName<<"', does the file exist?"<<std::endl;
0035       return false;
0036     }
0037 
0038     //Retrieve info, start with setup
0039     TsetupIn = (TTree*) RootInput->Get("Setup");
0040     if(!TsetupIn){
0041       std::cout<<"Could not retrieve the Setup tree, leaving"<<std::endl;
0042       return false;
0043     }
0044     setup=Setup::GetInstance();
0045     std::cout<<"Setup add "<<setup<<std::endl;
0046     //matchingbranch=TsetupIn->SetBranchAddress("setup",&setup);
0047     matchingbranch=TsetupIn->SetBranchAddress("setup",&rswptr);
0048     if(matchingbranch<0){
0049       std::cout<<"Error retrieving Setup info from the tree"<<std::endl;
0050       return false;
0051     }
0052     std::cout<<"Entries "<<TsetupIn->GetEntries()<<std::endl;
0053     TsetupIn->GetEntry(0);
0054     setup->Initialize(*rswptr);
0055     std::cout<<"Reading "<<RootInput->GetName()<<std::endl;
0056     std::cout<<"Setup Info "<<setup->IsInit()<<"  and  "<<setup->GetCellID(0,0)<<std::endl;
0057     //std::cout<<"Setup add now "<<setup<<std::endl;
0058 
0059     //Retrieve info, existing data?
0060     TdataIn = (TTree*) RootInput->Get("Data");
0061     if(!TdataIn){
0062       std::cout<<"Could not retrieve the Data tree, leaving"<<std::endl;
0063       return false;
0064     }
0065     matchingbranch=TdataIn->SetBranchAddress("event",&eventptr);
0066     if(matchingbranch<0){
0067       std::cout<<"Error retrieving Event info from the tree"<<std::endl;
0068       return false;
0069     }
0070     
0071     TcalibIn = (TTree*) RootInput->Get("Calib");
0072     if(!TcalibIn){
0073       std::cout<<"Could not retrieve Calib tree, leaving"<<std::endl;
0074       //return false;
0075     }
0076     else {
0077       matchingbranch=TcalibIn->SetBranchAddress("calib",&calibptr);
0078       if(matchingbranch<0){
0079         std::cout<<"Error retrieving calibration info from the tree"<<std::endl;
0080         TcalibIn=nullptr;
0081       }
0082     }
0083   } else {
0084     std::cout<<"An input file is required, aborting"<<std::endl;
0085     return false;
0086   }  
0087 
0088   return true;
0089 }
0090 
0091 // ****************************************************************************
0092 // Primary process function to start all calibrations
0093 // ****************************************************************************
0094 bool EventDisplay::Process(void){
0095   bool status;
0096   ROOT::EnableImplicitMT();
0097   
0098   // copy full calibration to different file and calculate energy
0099   status=Plot();
0100 
0101   return status;
0102 }
0103 
0104 
0105 //***********************************************************************************************
0106 //*********************** Calibration routine ***************************************************
0107 //***********************************************************************************************
0108 bool EventDisplay::Plot(){
0109   std::cout<<"Plotting events "<<plotEvt<<" to " << plotEvt+nEvts-1<<std::endl;
0110   if(plotMuonEvts) std::cout<<"Plotting muon-triggered events"<<std::endl;
0111 
0112   std::map<int,RunInfo> ri=readRunInfosFromFile(RunListInputName.Data(),debug,0);
0113   TdataIn->GetEntry(0);
0114   Int_t runNr = event.GetRunNumber();
0115   // Get run info object
0116   std::map<int,RunInfo>::iterator it=ri.find(runNr);
0117 
0118   int species = -1;
0119   species = GetSpeciesIntFromRunInfo(it->second);
0120   if (species == -1){
0121       std::cout << "WARNING: species unknown: " << it->second.species.Data() << "  aborting"<< std::endl;
0122       return false;
0123   }
0124   
0125   std::cout << "debug level set to : " << debug << std::endl;
0126   
0127   // create 3D histo
0128   // int towersx         = setup->GetNMaxColumn()+1;/*GetNMAxColumn() returns */
0129   float towersx_min   = setup->GetMinX();
0130   float towersx_max   = setup->GetMaxX();
0131   int towersx         = (towersx_max-towersx_min)/setup->GetCellWidth();
0132   float towersy_min   = setup->GetMinY();
0133   float towersy_max   = setup->GetMaxY();
0134   int towersy         = (towersy_max-towersy_min)/setup->GetCellHeight();
0135   float towersz_min   = setup->GetMinZ();
0136   float towersz_max   = setup->GetMaxZ();
0137   int towersz         = (towersz_max-towersz_min)/setup->GetCellDepth();
0138 
0139   TH3F*   hXYZMapEvt          = new TH3F("hXYZMapEvt","",towersz, towersz_min, towersz_max, towersx, towersx_min, towersx_max, towersy, towersy_min, towersy_max);
0140   TH3F*   hXYZMapEvt_Muon     = new TH3F("hXYZMapEvt_Muon","",towersz, towersz_min, towersz_max, towersx, towersx_min, towersx_max, towersy, towersy_min, towersy_max);
0141   TH3F*   hXYZMapEvt_nonMuon  = new TH3F("hXYZMapEvt_nonMuon","",towersz, towersz_min, towersz_max, towersx, towersx_min, towersx_max, towersy, towersy_min, towersy_max);
0142 
0143   TH1D* hX_energy_Evt    = new TH1D("hXenergyEvt","",towersx, towersx_min, towersx_max);
0144   hX_energy_Evt->Sumw2();
0145   TH1D* hY_energy_Evt    = new TH1D("hYenergyEvt","",towersy, towersy_min, towersy_max);
0146   hY_energy_Evt->Sumw2();
0147   TH1D* hZ_energy_Evt    = new TH1D("hZenergyEvt","",towersz, towersz_min, towersz_max);
0148   hZ_energy_Evt->Sumw2();
0149 
0150   TH1D* hX_energy_Evt_Muon    = new TH1D("hXenergyEvt_muon","",towersx, towersx_min, towersx_max);
0151   hX_energy_Evt_Muon->Sumw2();
0152   TH1D* hY_energy_Evt_Muon    = new TH1D("hYenergyEvt_muon","",towersy, towersy_min, towersy_max);
0153   hY_energy_Evt_Muon->Sumw2();
0154   TH1D* hZ_energy_Evt_Muon    = new TH1D("hZenergyEvt_muon","",towersz, towersz_min, towersz_max);
0155   hZ_energy_Evt_Muon->Sumw2();
0156 
0157   TH1D* hX_energy_Evt_nonMuon    = new TH1D("hXenergyEvt_nonMuon","",towersx, towersx_min, towersx_max);
0158   hX_energy_Evt_nonMuon->Sumw2();
0159   TH1D* hY_energy_Evt_nonMuon    = new TH1D("hYenergyEvt_nonMuon","",towersy, towersy_min, towersy_max);
0160   hY_energy_Evt_nonMuon->Sumw2();
0161   TH1D* hZ_energy_Evt_nonMuon    = new TH1D("hZenergyEvt_nonMuon","",towersz, towersz_min, towersz_max);
0162   hZ_energy_Evt_nonMuon->Sumw2();  
0163 
0164   // creating plotting directory
0165   StyleSettingsBasics(plotSuffix);
0166   TString outputDirPlots = GetPlotOutputDir();
0167   gSystem->Exec("mkdir -p "+outputDirPlots);
0168   
0169   TCanvas* canvas3D = new TCanvas("canvas3D","",0,0,1400,750);  // gives the page size
0170   DefaultCancasSettings( canvas3D, 0.12, 0.08, 0.05, 0.1);
0171   int evts=TdataIn->GetEntries();
0172   int evtsMuon= 0;
0173   for(int i=plotEvt; i<plotEvt+nEvts; i++){
0174     if(i > evts ) {
0175       std::cout<<"Requested event outside of the range. Aborting..."<<std::endl;
0176       return false;
0177     }
0178     TdataIn->GetEntry(i);
0179     if (i%100 == 0 && debug > 0) std::cout << "Reading " <<  i << " / " << evts << " events" << std::endl;
0180     double Etot = 0;
0181     int nCells  = 0;
0182     bool muontrigg = false;
0183     int locMuon = 0;
0184     std::map<int,Layer> layers;
0185     std::map<int, Layer>::iterator ithLayer;
0186     
0187     for(int j=0; j<event.GetNTiles(); j++){
0188       Caen* aTile=(Caen*)event.GetTile(j);
0189       if(aTile->GetE()>0.3 ){ 
0190         nCells++;
0191         int currLayer = setup->GetLayer(aTile->GetCellID());
0192         ithLayer=layers.find(currLayer);
0193         if(ithLayer!=layers.end()){
0194           ithLayer->second.nCells+=1;
0195           ithLayer->second.energy+=aTile->GetE();
0196           ithLayer->second.avX+=setup->GetX(aTile->GetCellID())*aTile->GetE();
0197           ithLayer->second.avY+=setup->GetY(aTile->GetCellID())*aTile->GetE();
0198         } else {
0199           layers[currLayer]=Layer();
0200           layers[currLayer].nCells+=1;
0201           layers[currLayer].energy+=aTile->GetE();
0202           layers[currLayer].avX+=setup->GetX(aTile->GetCellID())*aTile->GetE();
0203           layers[currLayer].avY+=setup->GetY(aTile->GetCellID())*aTile->GetE();
0204         }
0205       }
0206 
0207       if (aTile->GetLocalTriggerBit() == 1){
0208         locMuon++;        
0209       }      
0210     }
0211     
0212     if (nCells > 0) {
0213       int nLayerSingleCell = 0;
0214       for(ithLayer=layers.begin(); ithLayer!=layers.end(); ++ithLayer){
0215         if (ithLayer->second.nCells == 1) 
0216             nLayerSingleCell++;
0217       }
0218       double fracLayer1Cell = (double)nLayerSingleCell/(int)layers.size();
0219       double fracLocMuon = (double)locMuon/nCells;
0220       if (fracLocMuon > 0.6){ 
0221         muontrigg = true;
0222         evtsMuon++;
0223       }
0224       if (muontrigg && debug > 1){
0225           std::cout << "Event: " << i << "\tMuon triggered:\t" <<  fracLayer1Cell << "\t" << fracLocMuon << std::endl;
0226       }
0227       Float_t minE = 1e10;
0228       Float_t maxE = 0;
0229       for(int j=0; j<event.GetNTiles(); j++){
0230         Caen* aTile=(Caen*)event.GetTile(j);
0231         // remove bad channels from output
0232         Float_t energy = (Float_t)aTile->GetE();
0233         Etot+=aTile->GetE();
0234         if (energy < minE) minE = energy;
0235         if (energy > maxE) maxE = energy;
0236         if(energy>0.3){ 
0237           hXYZMapEvt->Fill(aTile->GetZ(),aTile->GetX(),aTile->GetY(),energy);
0238           hX_energy_Evt->Fill(aTile->GetX(), energy);
0239           hY_energy_Evt->Fill(aTile->GetY(), energy);
0240           hZ_energy_Evt->Fill(aTile->GetZ(), energy);
0241           if(debug > 2) std::cout << "Event, x, y, z, E: " << i << "\t" <<aTile->GetX()<< "\t" <<aTile->GetY()<< "\t" << aTile->GetZ()<< "\t" <<energy<<std::endl;
0242           if (aTile->GetLocalTriggerBit() == 1){
0243             hXYZMapEvt_Muon->Fill(aTile->GetZ(),aTile->GetX(),aTile->GetY(),energy);
0244             hX_energy_Evt_Muon->Fill(aTile->GetX(), energy);
0245             hY_energy_Evt_Muon->Fill(aTile->GetY(), energy);
0246             hZ_energy_Evt_Muon->Fill(aTile->GetZ(), energy);
0247 
0248           } else {
0249             hXYZMapEvt_nonMuon->Fill(aTile->GetZ(),aTile->GetX(),aTile->GetY(),energy);
0250             hX_energy_Evt_nonMuon->Fill(aTile->GetX(), energy);
0251             hY_energy_Evt_nonMuon->Fill(aTile->GetY(), energy);
0252             hZ_energy_Evt_nonMuon->Fill(aTile->GetZ(), energy);
0253           }
0254         } 
0255       }
0256 
0257       //**********************************************************************
0258       //********************* Plotting ***************************************
0259       //**********************************************************************  
0260       Float_t maxEX = 0;
0261       Float_t maxEY = 0;
0262       Float_t maxEZ = 0;
0263       for (Int_t k = 1; k < hX_energy_Evt->GetNbinsX()+1; k++){
0264         if (maxEX < hX_energy_Evt->GetBinContent(k)) maxEX = hX_energy_Evt->GetBinContent(k);
0265       }
0266       for (Int_t k = 1; k < hY_energy_Evt->GetNbinsX()+1; k++){
0267         if (maxEY < hY_energy_Evt->GetBinContent(k)) maxEY = hY_energy_Evt->GetBinContent(k);
0268       }
0269       for (Int_t k = 1; k < hZ_energy_Evt->GetNbinsX()+1; k++){
0270         if (maxEZ < hZ_energy_Evt->GetBinContent(k)) maxEZ = hZ_energy_Evt->GetBinContent(k);
0271       }
0272       //**********************************************************************
0273       // Create canvases for channel overview plotting
0274       //**********************************************************************
0275       Double_t textSizeRel = 0.035;
0276       StyleSettingsBasics("pdf");
0277       SetPlotStyle();
0278       if( (muontrigg&&plotMuonEvts) || !plotMuonEvts){
0279         EventDisplayWithSliceHighlighted( hXYZMapEvt, hX_energy_Evt, hY_energy_Evt, hZ_energy_Evt, 
0280                                          hXYZMapEvt_Muon, hX_energy_Evt_Muon, hY_energy_Evt_Muon, hX_energy_Evt_Muon, 
0281                                          hXYZMapEvt_nonMuon, hX_energy_Evt_nonMuon, hY_energy_Evt_nonMuon, hZ_energy_Evt_nonMuon, 
0282                                          i, Etot, maxE, maxEX, maxEY, maxEZ,  muontrigg,
0283                                          it->second, Form("%s/EventDisplay_muonHighlighed_evt", outputDirPlots.Data()), plotSuffix);    
0284       }
0285 
0286       if( (muontrigg&&plotMuonEvts) || !plotMuonEvts){
0287         EventDisplayWithSlice(  hXYZMapEvt, hX_energy_Evt, hY_energy_Evt, hZ_energy_Evt, 
0288                                 i, Etot, maxE, maxEX, maxEY, maxEZ,  muontrigg,
0289                                 it->second, Form("%s/EventDisplay_MonoChrome_evt", outputDirPlots.Data()), plotSuffix);    
0290 
0291         canvas3D->cd();
0292 
0293         SetStyleHistoTH3ForGraphs(hXYZMapEvt, "z", "x","y", 0.85*textSizeRel,textSizeRel, 0.85*textSizeRel,textSizeRel, 0.85*textSizeRel,textSizeRel, 1.1, 1.1, 1.15, 505, 510,510);
0294         hXYZMapEvt->SetMaximum(maxE);
0295         hXYZMapEvt->DrawCopy("box2z");
0296         DrawLatex(0.05, 0.94, GetStringFromRunInfo(it->second, 1), false, 0.85*textSizeRel, 42);
0297         if(muontrigg) DrawLatex(0.05, 0.90, Form("Event %d, muon triggered",i), false, 0.85*textSizeRel, 42);
0298         else DrawLatex(0.05, 0.90, Form("Event %d",i), false, 0.85*textSizeRel, 42);
0299         
0300         canvas3D->SaveAs( Form("%s/EventDisplay_Colored_evt%06i.%s", outputDirPlots.Data(), i, plotSuffix.Data()));
0301         canvas3D->ResetDrawn();
0302       }
0303     }
0304 
0305     hXYZMapEvt->Reset();
0306     hXYZMapEvt_Muon->Reset();
0307     hXYZMapEvt_nonMuon->Reset();
0308     hX_energy_Evt->Reset();
0309     hY_energy_Evt->Reset();
0310     hZ_energy_Evt->Reset();
0311     hX_energy_Evt_Muon->Reset();
0312     hY_energy_Evt_Muon->Reset();
0313     hZ_energy_Evt_Muon->Reset();
0314     hX_energy_Evt_nonMuon->Reset();
0315     hY_energy_Evt_nonMuon->Reset();
0316     hZ_energy_Evt_nonMuon->Reset();
0317   }
0318 
0319     delete hXYZMapEvt;
0320     delete hXYZMapEvt_Muon;
0321     delete hXYZMapEvt_nonMuon;
0322     delete hX_energy_Evt;
0323     delete hY_energy_Evt;
0324     delete hZ_energy_Evt;
0325     delete hX_energy_Evt_Muon;
0326     delete hY_energy_Evt_Muon;
0327     delete hZ_energy_Evt_Muon;
0328     delete hX_energy_Evt_nonMuon;
0329     delete hY_energy_Evt_nonMuon;
0330     delete hZ_energy_Evt_nonMuon;
0331     delete canvas3D;
0332     
0333     return true;
0334 }
[ 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!