Back to home page

EIC code displayed by LXR
 
 

    

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

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 "PlotHelper.h"
0017 #include "TRandom3.h"
0018 #include "TMath.h"
0019 #include "Math/Minimizer.h"
0020 #include "Math/MinimizerOptions.h"
0021 
0022 // ****************************************************************************
0023 // Checking and opening input and output files
0024 // ****************************************************************************
0025 bool 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       } else {
0082         std::cout<<"Loaded the calib object"<<std::endl;
0083       }
0084     }
0085   } else {
0086     std::cout<<"An input file is required, aborting"<<std::endl;
0087     return false;
0088   }  
0089 
0090   return true;
0091 }
0092 
0093 // ****************************************************************************
0094 // Primary process function to start all calibrations
0095 // ****************************************************************************
0096 bool EventDisplay::Process(void){
0097   bool status;
0098   ROOT::EnableImplicitMT();
0099   
0100   // copy full calibration to different file and calculate energy
0101   status=Plot();
0102 
0103   return status;
0104 }
0105 
0106 
0107 //***********************************************************************************************
0108 //*********************** Plotting routine ***************************************************
0109 //***********************************************************************************************
0110 bool EventDisplay::Plot(){
0111   std::cout<<"Plotting events "<<plotEvt<<" to " << plotEvt+nEvts-1<<std::endl;
0112   if(plotMuonEvts) std::cout<<"Plotting muon-triggered events"<<std::endl;
0113 
0114   std::map<int,RunInfo> ri=readRunInfosFromFile(RunListInputName.Data(),debug,0);
0115   TdataIn->GetEntry(0);
0116   Int_t runNr = event.GetRunNumber();
0117   // Get run info object
0118   std::map<int,RunInfo>::iterator it=ri.find(runNr);
0119 
0120   int species = -1;
0121   species = GetSpeciesIntFromRunInfo(it->second);
0122   if (species == -1){
0123       std::cout << "WARNING: species unknown: " << it->second.species.Data() << "  aborting"<< std::endl;
0124       return false;
0125   }
0126   
0127   std::cout << "debug level set to : " << debug << std::endl;
0128   
0129   // create 3D histo
0130   // int towersx         = setup->GetNMaxColumn()+1;/*GetNMAxColumn() returns */
0131   float towersx_min   = setup->GetMinX();
0132   float towersx_max   = setup->GetMaxX();
0133   int towersx         = (towersx_max-towersx_min)/setup->GetCellWidth();
0134   float towersy_min   = setup->GetMinY();
0135   float towersy_max   = setup->GetMaxY();
0136   int towersy         = (towersy_max-towersy_min)/setup->GetCellHeight();
0137   float towersz_min   = setup->GetMinZ();
0138   float towersz_max   = setup->GetMaxZ();
0139   int towersz         = (towersz_max-towersz_min)/setup->GetCellDepth();
0140 
0141   TH3F*   hXYZMapEvt          = new TH3F("hXYZMapEvt","",towersz, towersz_min, towersz_max, towersx, towersx_min, towersx_max, towersy, towersy_min, towersy_max);
0142   TH3F*   hXYZMapEvt_Muon     = new TH3F("hXYZMapEvt_Muon","",towersz, towersz_min, towersz_max, towersx, towersx_min, towersx_max, towersy, towersy_min, towersy_max);
0143   TH3F*   hXYZMapEvt_nonMuon  = new TH3F("hXYZMapEvt_nonMuon","",towersz, towersz_min, towersz_max, towersx, towersx_min, towersx_max, towersy, towersy_min, towersy_max);
0144 
0145   TH1D* hX_energy_Evt    = new TH1D("hXenergyEvt","",towersx, towersx_min, towersx_max);
0146   hX_energy_Evt->Sumw2();
0147   TH1D* hY_energy_Evt    = new TH1D("hYenergyEvt","",towersy, towersy_min, towersy_max);
0148   hY_energy_Evt->Sumw2();
0149   TH1D* hZ_energy_Evt    = new TH1D("hZenergyEvt","",towersz, towersz_min, towersz_max);
0150   hZ_energy_Evt->Sumw2();
0151 
0152   TH1D* hX_energy_Evt_Muon    = new TH1D("hXenergyEvt_muon","",towersx, towersx_min, towersx_max);
0153   hX_energy_Evt_Muon->Sumw2();
0154   TH1D* hY_energy_Evt_Muon    = new TH1D("hYenergyEvt_muon","",towersy, towersy_min, towersy_max);
0155   hY_energy_Evt_Muon->Sumw2();
0156   TH1D* hZ_energy_Evt_Muon    = new TH1D("hZenergyEvt_muon","",towersz, towersz_min, towersz_max);
0157   hZ_energy_Evt_Muon->Sumw2();
0158 
0159   TH1D* hX_energy_Evt_nonMuon    = new TH1D("hXenergyEvt_nonMuon","",towersx, towersx_min, towersx_max);
0160   hX_energy_Evt_nonMuon->Sumw2();
0161   TH1D* hY_energy_Evt_nonMuon    = new TH1D("hYenergyEvt_nonMuon","",towersy, towersy_min, towersy_max);
0162   hY_energy_Evt_nonMuon->Sumw2();
0163   TH1D* hZ_energy_Evt_nonMuon    = new TH1D("hZenergyEvt_nonMuon","",towersz, towersz_min, towersz_max);
0164   hZ_energy_Evt_nonMuon->Sumw2();  
0165 
0166   // creating plotting directory
0167   StyleSettingsBasics(plotSuffix);
0168   TString outputDirPlots = GetPlotOutputDir();
0169   gSystem->Exec("mkdir -p "+outputDirPlots);
0170   
0171   TCanvas* canvas3D = new TCanvas("canvas3D","",0,0,1400,750);  // gives the page size
0172   DefaultCancasSettings( canvas3D, 0.12, 0.08, 0.05, 0.1);
0173 
0174   // processing events
0175   int evts=TdataIn->GetEntries();
0176   int evtsMuon= 0;
0177   for(int i=plotEvt; i<plotEvt+nEvts; i++){
0178     if(i > evts ) {
0179       std::cout<<"Requested event outside of the range. Aborting..."<<std::endl;
0180       return false;
0181     }
0182     TdataIn->GetEntry(i);
0183     if (i%100 == 0 && debug > 0) std::cout << "Reading " <<  i << " / " << evts << " events" << std::endl;
0184     double Etot = 0;
0185     int nCells  = 0;
0186     bool muontrigg = false;
0187     int locMuon = 0;
0188     std::map<int,Layer> layers;
0189     std::map<int, Layer>::iterator ithLayer;
0190     
0191     for(int j=0; j<event.GetNTiles(); j++){
0192       if(plotHGCROC){
0193         Hgcroc* aTile = (Hgcroc*)event.GetTile(j);
0194         // no "energy" threshold for now - just plot everything
0195         nCells++;
0196         int currLayer   = setup->GetLayer(aTile->GetCellID());
0197         ithLayer        = layers.find(currLayer);
0198         if(ithLayer!=layers.end()){
0199           ithLayer->second.nCells+=1;
0200           double energy = 0;
0201           switch( dataTypeHGCROC ){
0202             case 0:
0203               energy = aTile->GetMaxSampleADC(); 
0204               if( debug>1 ) std::cout << "CellID: " << aTile->GetCellID() << "\t ADC: " << aTile->GetMaxSampleADC() << "\t energy " << energy << std::endl;
0205               break;
0206             case 1:
0207               energy = aTile->GetMaxSampleADC() - aTile->GetPedestal(); 
0208               if( debug>1 ) std::cout << "CellID: " << aTile->GetCellID() << "\t ADC: " << aTile->GetMaxSampleADC() << "\t pedestal (calib) " << calib.GetPedestalMeanH( aTile->GetCellID() ) << "\t pedestal (aTile) " << aTile->GetPedestal() << "\t energy " << energy << std::endl;
0209               break;
0210             case 2:
0211               energy = aTile->GetTOT(); 
0212               if( debug>1 ) std::cout << "CellID: " << aTile->GetCellID() << "\t TOT: " << aTile->GetTOT() << "\t energy " << energy << std::endl;
0213               break;
0214           }
0215           ithLayer->second.energy+=energy;
0216           ithLayer->second.avX+=setup->GetX(aTile->GetCellID())*energy;
0217           ithLayer->second.avY+=setup->GetY(aTile->GetCellID())*energy;
0218         } else {
0219           layers[currLayer]=Layer();
0220           layers[currLayer].nCells+=1;
0221           double energy = 0;
0222           switch( dataTypeHGCROC ){
0223             case 0:
0224               energy = aTile->GetMaxSampleADC(); 
0225               if( debug>1 ) std::cout << "CellID: " << aTile->GetCellID() << "\t ADC: " << aTile->GetMaxSampleADC() << "\t energy " << energy << std::endl;
0226               break;
0227             case 1:
0228               energy = aTile->GetMaxSampleADC() - aTile->GetPedestal(); 
0229               if( debug>1 ) std::cout << "CellID: " << aTile->GetCellID() << "\t ADC: " << aTile->GetMaxSampleADC() << "\t pedestal (calib) " << calib.GetPedestalMeanH( aTile->GetCellID() ) << "\t pedestal (aTile) " << aTile->GetPedestal() << "\t energy " << energy << std::endl;
0230               break;
0231             case 2:
0232               energy = aTile->GetTOT(); 
0233               if( debug>1 ) std::cout << "CellID: " << aTile->GetCellID() << "\t TOT: " << aTile->GetTOT() << "\t energy " << energy << std::endl;
0234               break;
0235           }
0236           layers[currLayer].energy+=energy;
0237           layers[currLayer].avX+=setup->GetX(aTile->GetCellID())*energy;
0238           layers[currLayer].avY+=setup->GetY(aTile->GetCellID())*energy;
0239         }
0240       } else {
0241         Caen* aTile=(Caen*)event.GetTile(j);
0242         if(aTile->GetE()>0.3 ){ 
0243           nCells++;
0244           int currLayer = setup->GetLayer(aTile->GetCellID());
0245           ithLayer=layers.find(currLayer);
0246           if(ithLayer!=layers.end()){
0247             ithLayer->second.nCells+=1;
0248             ithLayer->second.energy+=aTile->GetE();
0249             ithLayer->second.avX+=setup->GetX(aTile->GetCellID())*aTile->GetE();
0250             ithLayer->second.avY+=setup->GetY(aTile->GetCellID())*aTile->GetE();
0251           } else {
0252             layers[currLayer]=Layer();
0253             layers[currLayer].nCells+=1;
0254             layers[currLayer].energy+=aTile->GetE();
0255             layers[currLayer].avX+=setup->GetX(aTile->GetCellID())*aTile->GetE();
0256             layers[currLayer].avY+=setup->GetY(aTile->GetCellID())*aTile->GetE();
0257           }
0258         }
0259         if (aTile->GetLocalTriggerBit() == 1){
0260           locMuon++;        
0261         }      
0262       }
0263     }
0264     
0265     if (nCells > minTilesHit) {
0266       int nLayerSingleCell = 0;
0267       for(ithLayer=layers.begin(); ithLayer!=layers.end(); ++ithLayer){
0268         if (ithLayer->second.nCells == 1) 
0269             nLayerSingleCell++;
0270       }
0271       double fracLayer1Cell = (double)nLayerSingleCell/(int)layers.size();
0272       double fracLocMuon = (double)locMuon/nCells;
0273       if (fracLocMuon > 0.6){ 
0274         muontrigg = true;
0275         evtsMuon++;
0276       }
0277       if (muontrigg && debug > 1){
0278           std::cout << "Event: " << i << "\tMuon triggered:\t" <<  fracLayer1Cell << "\t" << fracLocMuon << std::endl;
0279       }
0280       Float_t minE = 1e10;
0281       Float_t maxE = 0;
0282       for(int j=0; j<event.GetNTiles(); j++){
0283         double energy  = 0;
0284         double x, y, z; 
0285         unsigned char muonTrigger = 0;
0286         if(plotHGCROC){
0287           Hgcroc* aTile   = (Hgcroc*)event.GetTile(j);
0288           if( calib.GetBadChannel(aTile->GetCellID())!=3) continue;
0289           switch( dataTypeHGCROC ){
0290             case 0:
0291               energy = aTile->GetMaxSampleADC();
0292               break;
0293             case 1:
0294               energy = aTile->GetMaxSampleADC()- calib.GetPedestalMeanL(aTile->GetCellID());
0295               break;
0296             case 2:
0297               energy = aTile->GetTOT();
0298               break;
0299           }
0300           x = (double) aTile->GetX();
0301           y = (double) aTile->GetY();
0302           z = (double) aTile->GetZ();
0303         } else {
0304           Caen* aTile=(Caen*)event.GetTile(j);
0305           // if( calib.GetBadChannel(aTile->GetCellID())!=3) continue; // commented out because of the calib object issues
0306           energy = (Float_t)aTile->GetE();
0307           x = (double) aTile->GetX();
0308           y = (double) aTile->GetY();
0309           z = (double) aTile->GetZ();
0310           muonTrigger = (unsigned char)aTile->GetLocalTriggerBit();
0311         }
0312         Etot+=energy;
0313         if (energy < minE) minE = energy;
0314         if (energy > maxE) maxE = energy;
0315         if(energy>0.3){ 
0316           hXYZMapEvt->Fill(z,x,y,energy);
0317           hX_energy_Evt->Fill(x, energy);
0318           hY_energy_Evt->Fill(y, energy);
0319           hZ_energy_Evt->Fill(z, energy);
0320           if(debug > 2) std::cout << "Event, x, y, z, E: " << i << "\t" <<x<< "\t" <<y<< "\t" << z<< "\t" <<energy<<std::endl;
0321           if (muonTrigger == 1){
0322             hXYZMapEvt_Muon->Fill(z,x,y,energy);
0323             hX_energy_Evt_Muon->Fill(x, energy);
0324             hY_energy_Evt_Muon->Fill(y, energy);
0325             hZ_energy_Evt_Muon->Fill(z, energy);
0326 
0327           } else {
0328             hXYZMapEvt_nonMuon->Fill(z,x,y,energy);
0329             hX_energy_Evt_nonMuon->Fill(x, energy);
0330             hY_energy_Evt_nonMuon->Fill(y, energy);
0331             hZ_energy_Evt_nonMuon->Fill(z, energy);
0332           }
0333         } 
0334       }
0335 
0336       //**********************************************************************
0337       //********************* Plotting ***************************************
0338       //**********************************************************************  
0339       Float_t maxEX = 0;
0340       Float_t maxEY = 0;
0341       Float_t maxEZ = 0;
0342       for (Int_t k = 1; k < hX_energy_Evt->GetNbinsX()+1; k++){
0343         if (maxEX < hX_energy_Evt->GetBinContent(k)) maxEX = hX_energy_Evt->GetBinContent(k);
0344       }
0345       for (Int_t k = 1; k < hY_energy_Evt->GetNbinsX()+1; k++){
0346         if (maxEY < hY_energy_Evt->GetBinContent(k)) maxEY = hY_energy_Evt->GetBinContent(k);
0347       }
0348       for (Int_t k = 1; k < hZ_energy_Evt->GetNbinsX()+1; k++){
0349         if (maxEZ < hZ_energy_Evt->GetBinContent(k)) maxEZ = hZ_energy_Evt->GetBinContent(k);
0350       }
0351       //**********************************************************************
0352       // Create canvases for channel overview plotting
0353       //**********************************************************************
0354       Double_t textSizeRel = 0.035;
0355       StyleSettingsBasics("pdf");
0356       SetPlotStyle();
0357       if( (muontrigg&&plotMuonEvts) || !plotMuonEvts){
0358         EventDisplayWithSliceHighlighted( hXYZMapEvt, hX_energy_Evt, hY_energy_Evt, hZ_energy_Evt, 
0359                                          hXYZMapEvt_Muon, hX_energy_Evt_Muon, hY_energy_Evt_Muon, hX_energy_Evt_Muon, 
0360                                          hXYZMapEvt_nonMuon, hX_energy_Evt_nonMuon, hY_energy_Evt_nonMuon, hZ_energy_Evt_nonMuon, 
0361                                          i, Etot, maxE, maxEX, maxEY, maxEZ,  muontrigg,
0362                                          it->second, Form("%s/EventDisplay_muonHighlighed_evt", outputDirPlots.Data()), plotSuffix);    
0363       }
0364 
0365       if( (muontrigg&&plotMuonEvts) || !plotMuonEvts){
0366         EventDisplayWithSlice(  hXYZMapEvt, hX_energy_Evt, hY_energy_Evt, hZ_energy_Evt, 
0367                                 i, Etot, maxE, maxEX, maxEY, maxEZ,  muontrigg,
0368                                 it->second, Form("%s/EventDisplay_MonoChrome_evt", outputDirPlots.Data()), plotSuffix);    
0369 
0370         canvas3D->cd();
0371 
0372         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);
0373         hXYZMapEvt->SetMaximum(maxE);
0374         hXYZMapEvt->DrawCopy("box2z");
0375         DrawLatex(0.05, 0.94, GetStringFromRunInfo(it->second, 1), false, 0.85*textSizeRel, 42);
0376         if(muontrigg) DrawLatex(0.05, 0.90, Form("Event %d, muon triggered",i), false, 0.85*textSizeRel, 42);
0377         else DrawLatex(0.05, 0.90, Form("Event %d",i), false, 0.85*textSizeRel, 42);
0378         
0379         canvas3D->SaveAs( Form("%s/EventDisplay_Colored_evt%06i.%s", outputDirPlots.Data(), i, plotSuffix.Data()));
0380         canvas3D->ResetDrawn();
0381       }
0382     }
0383 
0384     hXYZMapEvt->Reset();
0385     hXYZMapEvt_Muon->Reset();
0386     hXYZMapEvt_nonMuon->Reset();
0387     hX_energy_Evt->Reset();
0388     hY_energy_Evt->Reset();
0389     hZ_energy_Evt->Reset();
0390     hX_energy_Evt_Muon->Reset();
0391     hY_energy_Evt_Muon->Reset();
0392     hZ_energy_Evt_Muon->Reset();
0393     hX_energy_Evt_nonMuon->Reset();
0394     hY_energy_Evt_nonMuon->Reset();
0395     hZ_energy_Evt_nonMuon->Reset();
0396   }
0397 
0398   delete hXYZMapEvt;
0399   delete hXYZMapEvt_Muon;
0400   delete hXYZMapEvt_nonMuon;
0401   delete hX_energy_Evt;
0402   delete hY_energy_Evt;
0403   delete hZ_energy_Evt;
0404   delete hX_energy_Evt_Muon;
0405   delete hY_energy_Evt_Muon;
0406   delete hZ_energy_Evt_Muon;
0407   delete hX_energy_Evt_nonMuon;
0408   delete hY_energy_Evt_nonMuon;
0409   delete hZ_energy_Evt_nonMuon;
0410   delete canvas3D;
0411   
0412   return true;
0413 }
[ 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!