Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2026-04-21 07:51:39 

0001 #include "HGCROC_Waveform_Analysis.h"
0002 #include <vector>
0003 #include "TROOT.h"
0004 #include <bitset>
0005 #ifdef __APPLE__
0006 #include <unistd.h>
0007 #endif
0008 #include "TF1.h"
0009 #include "TFitResult.h"
0010 #include "TFitResultPtr.h"
0011 #include "TH1D.h"
0012 #include "TH2D.h"
0013 #include "TProfile.h"
0014 #include "TChain.h"
0015 #include "CommonHelperFunctions.h"
0016 #include "TileSpectra.h"
0017 #include "PlotHelper.h"
0018 #include "MultiCanvas.h"
0019 #include "TRandom3.h"
0020 #include "TMath.h"
0021 #include "Math/Minimizer.h"
0022 #include "Math/MinimizerOptions.h"
0023 
0024 #include "waveform_fitting/waveform_fit_base.h"
0025 #include "waveform_fitting/crystal_ball_fit.h"
0026 #include "waveform_fitting/max_sample_fit.h"
0027 
0028 // ****************************************************************************
0029 // Checking and opening input and output files
0030 // ****************************************************************************
0031 bool HGCROC_Waveform_Analysis::CheckAndOpenIO(void){
0032   int matchingbranch;
0033 
0034   //Need to check first input to get the setup...I do not think it is necessary
0035   std::cout <<"=============================================================" << std::endl;
0036   std::cout<<"Input name set to: "<<RootInputName.Data() <<std::endl;
0037   std::cout<<"Output name set to: "<<RootOutputName.Data() <<std::endl;
0038   std::cout<<"Calib name set to: "<<RootCalibInputName.Data() <<std::endl;
0039   std::cout <<"=============================================================" << std::endl;
0040   if(!RootInputName.IsNull()){
0041     //File exist?
0042     RootInput=new TFile(RootInputName.Data(),"READ");
0043     if(RootInput->IsZombie()){
0044       std::cout<<"Error opening '"<<RootInputName<<"', does the file exist?"<<std::endl;
0045       return false;
0046     }
0047 
0048     //Retrieve info, start with setup
0049     TsetupIn = (TTree*) RootInput->Get("Setup");
0050     if(!TsetupIn){
0051       std::cout<<"Could not retrieve the Setup tree, leaving"<<std::endl;
0052       return false;
0053     }
0054     setup=Setup::GetInstance();
0055     std::cout<<"Setup add "<<setup<<std::endl;
0056     //matchingbranch=TsetupIn->SetBranchAddress("setup",&setup);
0057     matchingbranch=TsetupIn->SetBranchAddress("setup",&rswptr);
0058     if(matchingbranch<0){
0059       std::cout<<"Error retrieving Setup info from the tree"<<std::endl;
0060       return false;
0061     }
0062     std::cout<<"Entries "<<TsetupIn->GetEntries()<<std::endl;
0063     TsetupIn->GetEntry(0);
0064     setup->Initialize(*rswptr);
0065     std::cout<<"Reading "<<RootInput->GetName()<<std::endl;
0066     std::cout<<"Setup Info "<<setup->IsInit()<<"  and  "<<setup->GetCellID(0,0)<<std::endl;
0067     //std::cout<<"Setup add now "<<setup<<std::endl;
0068 
0069     //Retrieve info, existing data?
0070     TdataIn = (TTree*) RootInput->Get("Data");
0071     if(!TdataIn){
0072       std::cout<<"Could not retrieve the Data tree, leaving"<<std::endl;
0073       return false;
0074     }
0075     matchingbranch=TdataIn->SetBranchAddress("event",&eventptr);
0076     if(matchingbranch<0){
0077       std::cout<<"Error retrieving Event info from the tree"<<std::endl;
0078       return false;
0079     }
0080     
0081     TcalibIn = (TTree*) RootInput->Get("Calib");
0082     if(!TcalibIn){
0083       std::cout<<"Could not retrieve Calib tree, leaving"<<std::endl;
0084       //return false;
0085     }
0086     else {
0087       matchingbranch=TcalibIn->SetBranchAddress("calib",&calibptr);
0088       if(matchingbranch<0){
0089         std::cout<<"Error retrieving calibration info from the tree"<<std::endl;
0090         TcalibIn=nullptr;
0091       }
0092     }
0093   }
0094     //Setup Output files
0095   if(!RootOutputName.IsNull()){    
0096     TString temp = RootOutputName;
0097     temp         = temp.ReplaceAll(".root","_Hists.root");
0098     SetRootOutputHists(temp);
0099     
0100     if (!IsAnalyseWaveForm && !IsExtractTimeWalk && !IsInvCrossTalk ){
0101       bool sCOF = CreateOutputRootFile();
0102       if (!sCOF) return false;
0103       
0104       TsetupOut = new TTree("Setup","Setup");
0105       
0106       //TsetupOut->Branch("setup",&setup);
0107       TsetupOut->Branch("setup",&rsw);
0108       TdataOut = new TTree("Data","Data");
0109       TdataOut->Branch("event",&event);
0110       TcalibOut = new TTree("Calib","Calib");
0111       TcalibOut->Branch("calib",&calib);
0112     }
0113   }
0114   
0115   if(!RootCalibInputName.IsNull()){
0116     RootCalibInput=new TFile(RootCalibInputName.Data(),"READ");
0117     if(RootCalibInput->IsZombie()){
0118       std::cout<<"Error opening '"<<RootCalibInputName<<"', does the file exist?"<<std::endl;
0119       return false;
0120     }
0121     TcalibIn = nullptr;
0122     TcalibIn = (TTree*) RootCalibInput->Get("Calib");
0123     if(!TcalibIn){
0124       std::cout<<"Could not retrieve Calib tree, leaving"<<std::endl;
0125       return false;
0126     } else {
0127       std::cout<<"Retrieved calib object from external input!"<<std::endl;
0128     }
0129     matchingbranch=TcalibIn->SetBranchAddress("calib",&calibptr);
0130     if(matchingbranch<0){
0131       std::cout<<"Error retrieving calibration info from the tree"<<std::endl;
0132       return false;
0133     }else {
0134       std::cout<<"Correctly set branch for external calib input."<<std::endl;
0135     }
0136   }
0137 
0138   std::cout <<"=============================================================" << std::endl;
0139   std::cout <<" Basic setup complete" << std::endl;
0140   std::cout <<"=============================================================" << std::endl;
0141   return true;    
0142 }
0143 
0144 // ****************************************************************************
0145 // Primary process function to start all calibrations
0146 // ****************************************************************************
0147 bool HGCROC_Waveform_Analysis::Process(void){
0148   bool status = true;
0149   ROOT::EnableImplicitMT();
0150   
0151   if (IsAnalyseWaveForm){
0152     status=AnalyseWaveForm();
0153   }
0154 
0155   if (IsExtractTimeWalk){
0156     status=ExtractTimeWalk();
0157   }
0158 
0159   if (IsInvCrossTalk){
0160     status=InvestigateCrossTalk();
0161   }
0162   
0163   return status;
0164 }
0165 
0166 // ****************************************************************************
0167 // Analyse waveform
0168 // ****************************************************************************
0169 bool HGCROC_Waveform_Analysis::AnalyseWaveForm(void){
0170   std::cout<<"Analyse Waveform"<<std::endl;
0171   int evts=TdataIn->GetEntries();
0172 
0173   std::map<int,RunInfo> ri=readRunInfosFromFile(RunListInputName.Data(),debug,0);
0174   
0175   std::map<int,TileSpectra> hSpectra;
0176   std::map<int, TileSpectra>::iterator ithSpectra;
0177   std::map<int,TileSpectra> hSpectraTrigg;
0178   std::map<int, TileSpectra>::iterator ithSpectraTrigg;
0179   TcalibIn->GetEntry(0);
0180   // check whether calib should be overwritten based on external text file
0181   if (OverWriteCalib){
0182     calib.ReadCalibFromTextFile(ExternalCalibFile,debug);
0183   }  
0184   if (ExternalToACalibOffSetFile.CompareTo("") != 0 ){
0185     calib.ReadExternalToAOffsets(ExternalToACalibOffSetFile,debug);
0186   }
0187 
0188   // reading first event in tree to extract runnumber & RO-type
0189   TdataIn->GetEntry(0);
0190   Int_t runNr           = event.GetRunNumber();
0191   ReadOut::Type typeRO  = event.GetROtype();
0192   std::cout<< "original run numbers calib: "<<calib.GetRunNumber() << "\t" << calib.GetRunNumberPed() << "\t" << calib.GetRunNumberMip() << std::endl;
0193   calib.SetRunNumber(runNr);
0194   calib.SetBeginRunTime(event.GetBeginRunTimeAlt());
0195   std::cout<< "reset run numbers calib: "<< calib.GetRunNumber() << "\t" << calib.GetRunNumberPed() << "\t" << calib.GetRunNumberMip() << std::endl;
0196   
0197   if (typeRO != ReadOut::Type::Hgcroc){
0198     std::cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl;
0199     std::cout << "This routine is not meant to run on CAEN data! Please check you are loading the correct inputs! \n ABORTING!!!" << std::endl;
0200     std::cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl;
0201     return false; 
0202   }
0203 
0204   //==================================================================================
0205   // create additional output hist
0206   //==================================================================================
0207   CreateOutputRootHistFile();
0208   
0209   int maxChannelPerLayer             = (setup->GetNMaxColumn()+1)*(setup->GetNMaxRow()+1)*(setup->GetNMaxModule()+1);
0210   int maxChannelPerLayerSingleMod    = (setup->GetNMaxColumn()+1)*(setup->GetNMaxRow()+1);
0211   TH2D* hHighestADCAbovePedVsLayer   = new TH2D( "hHighestEAbovePedVsLayer","Highest ADC above PED; layer; brd channel; #Sigma(ADC) (arb. units) ",
0212                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
0213   hHighestADCAbovePedVsLayer->SetDirectory(0);
0214   hHighestADCAbovePedVsLayer->Sumw2();
0215 
0216 
0217   TH2D* hspectraTOTvsCellID      = new TH2D( "hspectraTOTvsCellID","TOT spectrums CellID; cell ID; TOT (arb. units) ",
0218                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5, 4096,0,4096);
0219   hspectraTOTvsCellID->SetDirectory(0);
0220   TH2D* hspectraTOAvsCellID      = new TH2D( "hspectraTOAvsCellID","TOA spectrums CellID; cell ID; TOA (arb. units) ",
0221                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5, 1024,0,1024);
0222   hspectraTOAvsCellID->SetDirectory(0);
0223 
0224   TH2D* hSampleTOAVsCellID       = new TH2D( "hSampleTOAvsCellID","#sample ToA; cell ID; #sample TOA",
0225                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5, 20,0,20);
0226   hSampleTOAVsCellID->SetDirectory(0);
0227 
0228   TH2D* hTOAPsVsCellID       = new TH2D( "hTOAPsVsCellID","ToA (ns); cell ID; ToA (ps)",
0229                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5, 500,-250,0);
0230   hTOAPsVsCellID->SetDirectory(0);
0231   TH2D* hTOACorrNsVsCellID       = new TH2D( "hTOACorrNsVsCellID","ToA (ns); cell ID; ToA (ps)",
0232                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5, 500,-250,0);
0233   hTOACorrNsVsCellID->SetDirectory(0);
0234 
0235   TH2D* hSampleMaxADCVsCellID       = new TH2D( "hSampleMaxADCvsCellID","#sample ToA; cell ID; #sample Max ADC",
0236                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5, 20,0,20);
0237   hSampleMaxADCVsCellID->SetDirectory(0);
0238   TH2D* hSampleDiffvsCellID       = new TH2D( "hSampleDiffvsCellID","#sample diff; cell ID; #sample ToA-#sample Max ADC",
0239                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5, 8,-0.5,7.5);
0240   hSampleDiffvsCellID->SetDirectory(0);
0241   
0242   TH1D* hSampleTOA               = new TH1D( "hSampleTOA","#sample ToA; #sample TOA",
0243                                               20,0,20);
0244   TH1D* hSampleMaxADC            = new TH1D( "hSampleMaxADC","#sample ToA; #sample maxADC",
0245                                               20,0,20);
0246   TH1D* hSampleAboveTh           = new TH1D( "hSampleAboveTh","#sample ToA; #sample above th",
0247                                               20,0,20);
0248   TH1D* hSampleDiff              = new TH1D( "hSampleDiff","#sample ToA-#sample Max ADC; #sample maxADC-#sampleMax ADC",
0249                                               8,-0.5,7.5);
0250   TH1D* hSampleDiffMin           = new TH1D( "hSampleDiffMin","#sample ToA-#sample above th; #sample maxADC-#sample  above th",
0251                                               8,-0.5,7.5);
0252   
0253   RootOutputHist->mkdir("IndividualCells");
0254   RootOutputHist->mkdir("IndividualCellsTrigg");
0255   RootOutputHist->cd("IndividualCells");
0256   
0257   TH2D* h2DToAvsADC[setup->GetNMaxROUnit()+1][2];
0258   TProfile* hToAvsADC[setup->GetNMaxROUnit()+1][2];
0259   TH2D* h2DWaveFormHalfAsic[setup->GetNMaxROUnit()+1][2];
0260   TProfile* hWaveFormHalfAsic[setup->GetNMaxROUnit()+1][2];
0261   for (Int_t ro = 0; ro < setup->GetNMaxROUnit()+1; ro++){
0262     for (int h = 0; h< 2; h++){
0263       hToAvsADC[ro][h]    = new TProfile(Form("ADCTOA_Asic_%d_Half_%d",ro,h),
0264                                          Form("ADC-TOA correlation %d %d; TOA (arb. units); ADC (arb. units)",ro,h),1024/4,0,1024, "");
0265       h2DToAvsADC[ro][h]  = new TH2D(Form("h2DADCTOA_Asic_%d_Half_%d",ro,h),
0266                                      Form("2D ADC vs TOA %d %d; TOA (arb. units); ADC (arb. units)",ro,h), 1024/4,0,1024,1124,-100,1024);
0267       
0268       hWaveFormHalfAsic[ro][h]    = new TProfile(Form("WaveForm_Asic_%d_Half_%d",ro,h),
0269                                          Form("ADC-TOA correlation %d %d;  t (ns) ; ADC (arb. units)",ro,h),2200,-50,500);
0270       h2DWaveFormHalfAsic[ro][h]  = new TH2D(Form("h2DWaveForm_Asic_%d_Half_%d",ro,h),
0271                                         Form("2D ADC vs TOA %d %d;  t (ns) ; ADC (arb. units)",ro,h), 
0272                                         2200,-50,500, 1034, -10, 1024);
0273       
0274     }
0275   }
0276   
0277   ROOT::Math::MinimizerOptions::SetDefaultMinimizer("Minuit2", "Migrad");  
0278   if (maxEvents == -1){
0279     maxEvents = evts;
0280   } else {
0281     std::cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl;
0282     std::cout << "ATTENTION: YOU ARE RESETTING THE MAXIMUM NUMBER OF EVENTS TO BE PROCESSED TO: " << maxEvents << ". THIS SHOULD ONLY BE USED FOR TESTING!" << std::endl;
0283     std::cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl;
0284   }
0285   
0286   //==================================================================================
0287   // setup waveform builder for HGCROC data
0288   //==================================================================================
0289   waveform_fit_base *waveform_builder = nullptr;
0290   waveform_builder = new max_sample_fit();
0291   double minNSigma = 5;
0292   int nCellsAboveTOA  = 0;
0293   int nCellsAboveMADC = 0;
0294   double totADCs      = 0.;
0295   //==================================================================================
0296   // process events from tree
0297   //==================================================================================
0298   for(int i=0; i<evts && i < maxEvents ; i++){
0299     if (i%5000 == 0&& i > 0 && debug > 0) std::cout << "Reading " <<  i << "/" << evts << " events"<< std::endl;
0300     if (debug > 2 ){
0301       std::cout << "************************************* NEW EVENT " << i << "  *********************************" << std::endl;
0302     }
0303     TdataIn->GetEntry(i);
0304    
0305     // initialize counters per event
0306     nCellsAboveTOA  = 0;
0307     nCellsAboveMADC = 0;
0308     totADCs         = 0.;
0309     //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
0310     // process tiles in event
0311     //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++    
0312     for(int j=0; j<event.GetNTiles(); j++){
0313       //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
0314       // histo filling for CAEN specific things
0315       //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++      
0316       if (typeRO == ReadOut::Type::Caen) {
0317         // nothing to be done
0318         return false;
0319       //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++  
0320       // histo filling for HGCROC specific things & resetting of ped
0321       //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++        
0322       } else if (typeRO == ReadOut::Type::Hgcroc) {
0323         Hgcroc* aTile=(Hgcroc*)event.GetTile(j);
0324         int cellID = aTile->GetCellID();
0325         ithSpectra=hSpectra.find(cellID);
0326         ithSpectraTrigg=hSpectraTrigg.find(cellID);
0327         
0328         // get pedestal values from calib object
0329         double ped    = calib.GetPedestalMeanL(cellID);
0330         double pedErr = calib.GetPedestalSigL(cellID);
0331         if (ped == -1000){
0332           ped     = calib.GetPedestalMeanH(cellID);
0333           pedErr  = calib.GetPedestalSigH(cellID);
0334           if (ped == -1000){
0335             ped     = aTile->GetPedestal();
0336             pedErr  = 5;
0337           }
0338         }
0339         // reevaluate waveform
0340         waveform_builder->set_waveform(aTile->GetADCWaveform());
0341         waveform_builder->fit_with_average_ped(ped);
0342         aTile->SetIntegratedADC(waveform_builder->get_E());
0343         aTile->SetPedestal(waveform_builder->get_pedestal());
0344         
0345         // obtain integrated tile quantities for histo filling
0346         double adc      = aTile->GetIntegratedADC();
0347         double tot      = aTile->GetRawTOT();
0348         double toa      = aTile->GetRawTOA();
0349         bool hasTOA     = false;
0350         bool hasTOT     = false;
0351         
0352         // check if we have a toa
0353         if (toa > 0)
0354           hasTOA        = true;
0355         // check if we have a tot
0356         if (tot > 0)
0357           hasTOT        = true;
0358         Int_t nSampTOA  = aTile->GetFirstTOASample();        
0359         double toaNs      = (double)aTile->GetLinearizedRawTOA()*aTile->GetTOATimeResolution();
0360         double toaCorrNs  = toaNs;
0361         if (calib.GetToAOff(cellID) != -1000.){
0362           // correct ToA sample and return correct nSampleTOA
0363           nSampTOA      = aTile->SetCorrectedTOA(calib.GetToAOff(cellID)); 
0364           toaCorrNs     = (double)(aTile->GetCorrectedTOA())*aTile->GetTOATimeResolution();;
0365         }
0366         totADCs         = totADCs+adc;
0367         // obtain samples in which TOA, TOT, maximum ADC are firing
0368         
0369         Int_t nMaxADC   = aTile->GetMaxSampleADC();
0370         Int_t nADCFirst = aTile->GetFirstSampleAboveTh();
0371         Int_t nDiffFirstT= nSampTOA-nADCFirst;
0372         Int_t nDiffFirstM= nMaxADC-nADCFirst;
0373         Int_t nDiffMaxT  = nMaxADC-nSampTOA;
0374         Int_t nOffEmpty  = 2;
0375         
0376         if (adc > 10)
0377           nCellsAboveMADC++; 
0378         
0379         if (hasTOA){
0380           nCellsAboveTOA++;
0381           hSampleTOA->Fill(nSampTOA);
0382           hSampleTOAVsCellID->Fill(cellID,nSampTOA);
0383         }
0384         if (adc > ped+2*pedErr){
0385           hSampleMaxADC->Fill(nMaxADC);
0386           hSampleMaxADCVsCellID->Fill(cellID,nMaxADC);
0387           hSampleAboveTh->Fill(nADCFirst);
0388         }
0389         
0390         if (hasTOA && adc > ped+2*pedErr){
0391           hSampleDiff->Fill(nSampTOA-nMaxADC);
0392           hSampleDiffvsCellID->Fill(cellID,nSampTOA-nMaxADC);
0393           hSampleDiffMin->Fill(nSampTOA-nADCFirst);
0394         }
0395         
0396         if(hasTOA)  hspectraTOAvsCellID->Fill(cellID,toa);
0397         if(hasTOT)  hspectraTOTvsCellID->Fill(cellID,tot);
0398 
0399         Int_t asic      = setup->GetROunit(cellID);
0400         Int_t roCh      = setup->GetROchannel(cellID);
0401         bool fillToACorr = false;
0402         if(fixedTOASample > -1 && fixedTOASample == nSampTOA)
0403           fillToACorr = true;
0404         else if(fixedTOASample == -1 )
0405           fillToACorr = true;
0406         
0407         if (hasTOA){
0408           hTOAPsVsCellID->Fill(cellID,toaNs);
0409           hTOACorrNsVsCellID->Fill(cellID,toaCorrNs);
0410           // fill overview hists for half asics
0411           if( fillToACorr ){   
0412             hToAvsADC[asic][int(roCh/36)]->Fill(toa, adc);
0413             h2DToAvsADC[asic][int(roCh/36)]->Fill(toa, adc);
0414             for (int k = 0; k < (int)(aTile->GetADCWaveform()).size(); k++ ){
0415               double timetemp = ((k+nOffEmpty)*1024 + aTile->GetCorrectedTOA())*aTile->GetTOATimeResolution() ;
0416               hWaveFormHalfAsic[asic][int(roCh/36)]->Fill(timetemp,(aTile->GetADCWaveform()).at(k)-ped);
0417               h2DWaveFormHalfAsic[asic][int(roCh/36)]->Fill(timetemp,(aTile->GetADCWaveform()).at(k)-ped);
0418             }
0419           }
0420         }
0421         
0422         int layer     = setup->GetLayer(cellID);
0423         int chInLayer = setup->GetChannelInLayerFull(cellID);
0424         // Detailed debug info with print of waveform
0425         if (debug > 3){
0426           aTile->PrintWaveFormDebugInfo(calib.GetPedestalMeanH(cellID), calib.GetPedestalMeanL(cellID), calib.GetPedestalSigL(cellID));
0427         }
0428       
0429         // fill spectra histos
0430         if(ithSpectra!=hSpectra.end()){
0431           ithSpectra->second.FillExtHGCROC(adc,toa,tot,nSampTOA,fixedTOASample);
0432           if(fillToACorr) ithSpectra->second.FillWaveformVsTime(aTile->GetADCWaveform(), aTile->GetCorrectedTOA(), calib.GetPedestalMeanH(cellID),nOffEmpty);
0433         } else {
0434           RootOutputHist->cd("IndividualCells");
0435           hSpectra[cellID]=TileSpectra("full",3,cellID,calib.GetTileCalib(cellID),event.GetROtype(),debug);
0436           hSpectra[cellID].FillExtHGCROC(adc,toa,tot,nSampTOA, fixedTOASample);
0437           if(fillToACorr) hSpectra[cellID].FillWaveformVsTime(aTile->GetADCWaveform(), aTile->GetCorrectedTOA(), calib.GetPedestalMeanH(cellID),nOffEmpty);
0438         }
0439         // fill spectra histos for signals with ToA
0440         if (hasTOA){      // needed for summing tests
0441           hHighestADCAbovePedVsLayer->Fill(layer,chInLayer, adc);
0442           if(ithSpectraTrigg!=hSpectraTrigg.end()){
0443             ithSpectraTrigg->second.FillExtHGCROC(adc,toa,tot,nSampTOA,fixedTOASample);
0444             if(fillToACorr) ithSpectraTrigg->second.FillWaveformVsTime(aTile->GetADCWaveform(), aTile->GetCorrectedTOA(), calib.GetPedestalMeanH(cellID),nOffEmpty);
0445           } else {
0446             RootOutputHist->cd("IndividualCellsTrigg");
0447             hSpectraTrigg[cellID]=TileSpectra("signal",3,cellID,calib.GetTileCalib(cellID),event.GetROtype(),debug);
0448             hSpectraTrigg[cellID].FillExtHGCROC(adc,toa,tot,nSampTOA,fixedTOASample);
0449             if(fillToACorr) hSpectraTrigg[cellID].FillWaveformVsTime(aTile->GetADCWaveform(), aTile->GetCorrectedTOA(), calib.GetPedestalMeanH(cellID),nOffEmpty);
0450           }
0451         }
0452       }
0453     }
0454   }
0455 
0456   //==================================================================================
0457   // Setup general plotting infos
0458   //==================================================================================    
0459   std::map<int,RunInfo>::iterator it=ri.find(runNr);
0460   TString outputDirPlots = GetPlotOutputDir();
0461   Double_t textSizeRel = 0.035;
0462 
0463   gSystem->Exec("mkdir -p "+outputDirPlots);
0464   TString toaSampleAdd  = "";
0465   TString outputDirPlotsMore = outputDirPlots;
0466   if (fixedTOASample != -1){
0467     ExtPlot             = 1;
0468     toaSampleAdd        = Form("/nTOA_%02d",fixedTOASample);
0469     outputDirPlotsMore  = Form("%s%s", outputDirPlots.Data(),toaSampleAdd.Data());
0470     gSystem->Exec("mkdir -p "+outputDirPlotsMore);
0471   }
0472 
0473   
0474   StyleSettingsBasics("pdf");
0475   SetPlotStyle();  
0476   
0477   //==================================================================================
0478   // Create Summary plots
0479   //==================================================================================    
0480   TCanvas* canvas2DSigQA = new TCanvas("canvas2DSigQA","",0,0,1450,1300);  // gives the page size
0481   DefaultCanvasSettings( canvas2DSigQA, 0.08, 0.13, 0.045, 0.07);
0482 
0483   PlotSimple2D( canvas2DSigQA, hSampleTOAVsCellID, (double)it->second.samples,setup->GetMaxCellID()+1, textSizeRel, Form("%s/SampleTOAvsCellID.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
0484   PlotSimple2D( canvas2DSigQA, hSampleMaxADCVsCellID, (double)it->second.samples, setup->GetMaxCellID()+1, textSizeRel, Form("%s/SampleMaxADCvsCellID.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
0485   PlotSimple2D( canvas2DSigQA, hSampleDiffvsCellID, -10000, setup->GetMaxCellID()+1, textSizeRel, Form("%s/SampleDiffvsCellID.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
0486   PlotSimple2D( canvas2DSigQA, hspectraTOAvsCellID, -10000, setup->GetMaxCellID()+1, textSizeRel, Form("%s/TOAvsCellID.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
0487   PlotSimple2D( canvas2DSigQA, hspectraTOTvsCellID, -10000, setup->GetMaxCellID()+1, textSizeRel, Form("%s/TOTvsCellID.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
0488 
0489   PlotSimple2D( canvas2DSigQA, hTOAPsVsCellID, -10000, setup->GetMaxCellID()+1, textSizeRel, Form("%s/TOAPsvsCellID.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
0490   PlotSimple2D( canvas2DSigQA, hTOACorrNsVsCellID, -10000, setup->GetMaxCellID()+1, textSizeRel, Form("%s/TOACorrPsvsCellID.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
0491   
0492   canvas2DSigQA->SetLogz(1);
0493   PlotSimple2D( canvas2DSigQA, hHighestADCAbovePedVsLayer, -10000, -10000, textSizeRel, Form("%s/MaxADCAboveNoise_vsLayer.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);    
0494   
0495   Int_t diffBins = 4;
0496   for (Int_t ro = 0; ro < setup->GetNMaxROUnit()+1; ro++){
0497     for (int h = 0; h< 2; h++){
0498       Double_t largestDiff = 0.;
0499       Int_t bin         = 1;
0500       for (Int_t b = 1; b < hToAvsADC[ro][h]->GetNbinsX(); b++){
0501         Int_t maxBin  = (b+diffBins);
0502         if (maxBin > hToAvsADC[ro][h]->GetNbinsX()){
0503           maxBin      = maxBin-hToAvsADC[ro][h]->GetNbinsX();
0504         }
0505         Double_t diff = hToAvsADC[ro][h]->GetBinContent(maxBin) - hToAvsADC[ro][h]->GetBinContent(b);
0506         if (largestDiff < diff){
0507           largestDiff = diff;
0508           bin         = b;
0509         }
0510       }
0511       Int_t centerbin = bin+Int_t(diffBins/2);
0512       if (centerbin > (hToAvsADC[ro][h]->GetNbinsX()-diffBins/2))
0513         centerbin     = centerbin-(hToAvsADC[ro][h]->GetNbinsX()-diffBins/2);
0514       
0515       Int_t offset = (Int_t)(hToAvsADC[ro][h]->GetBinCenter(centerbin));
0516       std::cout << ro <<"\t" <<h << "\t" << offset << std::endl;
0517       
0518       Plot2DWithProfile( canvas2DSigQA, h2DToAvsADC[ro][h], hToAvsADC[ro][h], -10000, -10000, textSizeRel,
0519                         Form("%s/ToAvsADC_Asic_%d_Half_%d.%s", outputDirPlotsMore.Data(), ro, h, plotSuffix.Data()), 
0520                         it->second, 1, kFALSE, "colz",true, Form("Asic %d, Half %d",ro,h), offset);
0521       
0522       h2DWaveFormHalfAsic[ro][h]->GetXaxis()->SetRangeUser(-25, (it->second.samples)*25);
0523       Plot2DWithProfile( canvas2DSigQA, h2DWaveFormHalfAsic[ro][h], hWaveFormHalfAsic[ro][h], -10000, -10000, textSizeRel,
0524                         Form("%s/Waveform_Asic_%d_Half_%d.%s", outputDirPlotsMore.Data(), ro, h, plotSuffix.Data()), 
0525                         it->second, 1, kFALSE, "colz",true, Form("Asic %d, Half %d",ro,h), -1);
0526     }
0527   }
0528   
0529   TCanvas* canvas1DSimple = new TCanvas("canvas1DSimple","",0,0,1450,1300);  // gives the page size
0530   DefaultCanvasSettings( canvas1DSimple, 0.08, 0.03, 0.03, 0.07);
0531 
0532   PlotSimple1D(canvas1DSimple, hSampleTOA, -10000, (double)it->second.samples, textSizeRel, Form("%s/NSampleToA.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1);
0533   PlotSimple1D(canvas1DSimple, hSampleMaxADC, -10000, (double)it->second.samples, textSizeRel, Form("%s/NSampleMaxADC.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1);
0534   PlotSimple1D(canvas1DSimple, hSampleAboveTh, -10000, (double)it->second.samples, textSizeRel, Form("%s/NSampleAboveTh.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1);
0535   PlotSimple1D(canvas1DSimple, hSampleDiff, -10000, (double)it->second.samples, textSizeRel, Form("%s/NSampleDiff.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1);
0536   PlotSimple1D(canvas1DSimple, hSampleDiffMin, -10000, (double)it->second.samples, textSizeRel, Form("%s/NSampleDiffMin.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1);
0537 
0538   // print general calib info
0539   calib.PrintGlobalInfo();  
0540 
0541   //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++  
0542   // Plotting Single layers
0543   //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++  
0544   DetConf::Type detConf = setup->GetDetectorConfig();
0545   
0546   // set pixel text size
0547   Int_t textSizePixel = 30;
0548   MultiCanvas panelPlot(detConf, "AnaWave");
0549   bool init1D = panelPlot.Initialize(1);
0550   MultiCanvas panelPlot2D(detConf, "AnaWave2D");
0551   bool init2D = panelPlot2D.Initialize(2);
0552     
0553   panelPlot2D.PlotCorr2DLayer(hSpectra, 1, -25, (it->second.samples)*25, 0, 300,
0554                               Form("%s/Waveform",outputDirPlots.Data()), plotSuffix.Data(), it->second, &calib );
0555   panelPlot2D.PlotCorr2DLayer(hSpectra, 2, 0, 1024, 0, 300,
0556                               Form("%s/TOA_ADC",outputDirPlots.Data()), plotSuffix.Data(), it->second, &calib );
0557   panelPlot2D.PlotCorr2DLayer(hSpectra, 3, 0, 1024, 0, it->second.samples,
0558                               Form("%s/TOA_Sample",outputDirPlots.Data()), plotSuffix.Data(), it->second, &calib );
0559   panelPlot2D.PlotCorr2DLayer(hSpectraTrigg, 1, -25, (it->second.samples)*25, 0, 300,
0560                               Form("%s/WaveformSignal",outputDirPlots.Data()), plotSuffix.Data(), it->second, &calib );
0561   panelPlot2D.PlotCorr2DLayer(hSpectraTrigg, 2, 0, 1024, 0, 300,
0562                               Form("%s/TOA_ADC_Signal",outputDirPlots.Data()), plotSuffix.Data(), it->second, &calib );
0563   panelPlot2D.PlotCorr2DLayer(hSpectraTrigg, 3, 0, 1024, 0, it->second.samples,
0564                               Form("%s/TOA_Sample_Signal",outputDirPlots.Data()), plotSuffix.Data(), it->second, &calib );
0565   if (ExtPlot > 1){
0566     panelPlot.PlotSpectra( hSpectra, 0, -100, 1024, 1.2,
0567                            Form("%s/Spectra_ADC" ,outputDirPlots.Data()), plotSuffix.Data(), it->second, &calib);
0568     panelPlot.PlotSpectra( hSpectra, 3, 0, 1024, 1.2,
0569                            Form("%s/TOA" ,outputDirPlots.Data()), plotSuffix.Data(), it->second, &calib);
0570     panelPlot.PlotSpectra( hSpectra, 4, 0, 4096, 1.2,
0571                            Form("%s/TOT" ,outputDirPlots.Data()), plotSuffix.Data(), it->second, &calib);
0572   }
0573     
0574   //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++  
0575   // Saving histo putput
0576   //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++  
0577   std::cout << "Writing output hist file" << std::endl;
0578   RootOutputHist->cd();
0579   hHighestADCAbovePedVsLayer->Write();
0580   hSampleTOA->Write();
0581   hSampleMaxADC->Write();
0582   hSampleTOAVsCellID->Write();
0583   hSampleMaxADCVsCellID->Write();
0584   hspectraTOAvsCellID->Write();
0585   hspectraTOTvsCellID->Write();
0586 
0587   for (Int_t ro = 0; ro < setup->GetNMaxROUnit()+1; ro++){
0588     for (int h = 0; h< 2; h++){
0589       hToAvsADC[ro][h]->Write();
0590       h2DToAvsADC[ro][h]->Write();
0591     }
0592   }
0593   
0594   RootOutputHist->cd("IndividualCells");
0595   for(ithSpectra=hSpectra.begin(); ithSpectra!=hSpectra.end(); ++ithSpectra){
0596     ithSpectra->second.WriteExt(true);
0597   }
0598   RootOutputHist->cd("IndividualCellsTrigg");
0599   if (typeRO == ReadOut::Type::Hgcroc){
0600     for(ithSpectraTrigg=hSpectraTrigg.begin(); ithSpectraTrigg!=hSpectraTrigg.end(); ++ithSpectraTrigg){
0601       ithSpectraTrigg->second.WriteExt(true);
0602     }
0603   }
0604 
0605   RootInput->Close();
0606   return true;
0607 }
0608 
0609 // ****************************************************************************
0610 // Extract Time Walk
0611 // ****************************************************************************
0612 bool HGCROC_Waveform_Analysis::ExtractTimeWalk(void){
0613   std::cout<<"Analyse Waveform"<<std::endl;
0614   int evts=TdataIn->GetEntries();
0615 
0616   std::map<int,RunInfo> ri=readRunInfosFromFile(RunListInputName.Data(),debug,0);
0617   
0618   std::map<int,TileSpectra> hSpectra;
0619   std::map<int, TileSpectra>::iterator ithSpectra;
0620   std::map<int,TileSpectra> hSpectraAlt;
0621   std::map<int, TileSpectra>::iterator ithSpectraAlt;
0622   // set calib entry pointer
0623   TcalibIn->GetEntry(0);
0624 
0625   // obtain runnumber & RO-type
0626   TdataIn->GetEntry(0);
0627   Int_t runNr          = event.GetRunNumber();
0628   ReadOut::Type typeRO = event.GetROtype();
0629   std::map<int,RunInfo>::iterator it=ri.find(runNr);
0630   
0631   //==================================================================================
0632   // create additional output hist
0633   //==================================================================================
0634   CreateOutputRootHistFile();
0635 
0636   TH2D* hTOACorrNsVsCellID        = new TH2D( "hTOACorrNsVsCellID","ToA (ns); cell ID; ToA (ns)",
0637                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5, 500,-250,0);
0638   hTOACorrNsVsCellID->SetDirectory(0);
0639   TH2D* hTOACorrVsCellID          = new TH2D( "hTOACorrVsCellID"," ; cell ID; linearized ToA",
0640                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5, 1280,-10*1024,0);
0641   hTOACorrVsCellID->SetDirectory(0);
0642 
0643   TH1D* hSampleTOA               = new TH1D( "hSampleTOA","#sample ToA; #sample TOA",
0644                                               it->second.samples,0,it->second.samples);
0645   TH1D* hSampleMaxADC            = new TH1D( "hSampleMaxADC","#sample ToA; #sample maxADC",
0646                                               it->second.samples,0,it->second.samples);
0647   TH1D* hSampleAboveTh           = new TH1D( "hSampleAboveTh","#sample ToA; #sample above th",
0648                                               it->second.samples,0,it->second.samples);
0649   
0650   RootOutputHist->mkdir("IndividualCells");
0651   RootOutputHist->mkdir("IndividualCellsAlt");
0652   RootOutputHist->cd("IndividualCells");
0653   
0654   ROOT::Math::MinimizerOptions::SetDefaultMinimizer("Minuit2", "Migrad");  
0655   if (maxEvents == -1){
0656     maxEvents = evts;
0657   } else {
0658     std::cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl;
0659     std::cout << "ATTENTION: YOU ARE RESETTING THE MAXIMUM NUMBER OF EVENTS TO BE PROCESSED TO: " << maxEvents << ". THIS SHOULD ONLY BE USED FOR TESTING!" << std::endl;
0660     std::cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl;
0661   }
0662   
0663   //==================================================================================
0664   // process events from tree
0665   //==================================================================================
0666   for(int i=0; i<evts && i < maxEvents ; i++){
0667     if (i%5000 == 0&& i > 0 && debug > 0) std::cout << "Reading " <<  i << "/" << evts << " events"<< std::endl;
0668     if (debug > 2 && typeRO == ReadOut::Type::Hgcroc){
0669       std::cout << "************************************* NEW EVENT " << i << "  *********************************" << std::endl;
0670     }
0671     TdataIn->GetEntry(i);
0672     //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
0673     // process tiles in event
0674     //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++    
0675     for(int j=0; j<event.GetNTiles(); j++){
0676       //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++  
0677       // histo filling for HGCROC specific things & resetting of ped
0678       //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++        
0679       Hgcroc* aTile=(Hgcroc*)event.GetTile(j);
0680       int cellID = aTile->GetCellID();
0681       ithSpectra=hSpectra.find(cellID);
0682       ithSpectraAlt=hSpectraAlt.find(cellID);
0683       
0684       // get pedestal values from calib object
0685       double ped    = aTile->GetPedestal();
0686       
0687       // obtain integrated tile quantities for histo filling
0688       double adc      = aTile->GetIntegratedADC();
0689       double toaRaw   = aTile->GetRawTOA();
0690       double toa      = (double)(aTile->GetCorrectedTOA());
0691       Int_t nSampTOA  = aTile->GetCorrectedFirstTOASample(calib.GetToAOff(cellID));
0692       bool hasTOA     = false;
0693       
0694       // check if we have a toa
0695       if (toaRaw > 0)
0696         hasTOA        = true;
0697 
0698       double toaNs    = toa*aTile->GetTOATimeResolution();;
0699       
0700       Int_t nMaxADC   = aTile->GetMaxSampleADC();
0701       Int_t nADCFirst = aTile->GetFirstSampleAboveTh();
0702       Int_t nOffEmpty  = 2;
0703       
0704       hSampleMaxADC->Fill(nMaxADC);
0705       hSampleAboveTh->Fill(nADCFirst);
0706       
0707       if (hasTOA){
0708         hSampleTOA->Fill(nSampTOA);
0709         hTOACorrNsVsCellID->Fill(cellID,toaNs);
0710         hTOACorrVsCellID->Fill(cellID,toa);
0711       }
0712               
0713       int layer     = setup->GetLayer(cellID);
0714       int chInLayer = setup->GetChannelInLayerFull(cellID);
0715       // Detailed debug info with print of waveform
0716       if (debug > 3){
0717         aTile->PrintWaveFormDebugInfo(calib.GetPedestalMeanH(cellID), calib.GetPedestalMeanL(cellID), calib.GetPedestalSigL(cellID));
0718       }
0719     
0720       // fill spectra histos
0721       if(ithSpectra!=hSpectra.end()){
0722         ithSpectra->second.FillWaveformVsTime(aTile->GetADCWaveform(), aTile->GetCorrectedTOA(), ped, nOffEmpty);
0723       } else {
0724         RootOutputHist->cd("IndividualCells");
0725         hSpectra[cellID]=TileSpectra("full",6,cellID,calib.GetTileCalib(cellID),event.GetROtype(),debug);
0726         hSpectra[cellID].FillWaveformVsTime(aTile->GetADCWaveform(), aTile->GetCorrectedTOA(), ped, nOffEmpty);
0727       }
0728       // fill spectra histos for signals with ToA
0729       if (hasTOA){      // needed for summing tests
0730         Double_t tempADC    = 0.;
0731         Int_t nSampleTOARaw = aTile->GetFirstTOASample(); 
0732         if (nMaxADC > 0 && nMaxADC < (int)it->second.samples){
0733           tempADC = (aTile->GetADCWaveform()).at(nMaxADC);
0734         }
0735         if(ithSpectraAlt!=hSpectraAlt.end()){
0736           ithSpectraAlt->second.FillExtHGCROC(tempADC,toa,-1,nSampTOA,-1);
0737           ithSpectraAlt->second.FillMaxVsTime(tempADC, aTile->GetCorrectedTOA(),nOffEmpty, nMaxADC);
0738         } else {
0739           RootOutputHist->cd("IndividualCellsAlt");
0740           hSpectraAlt[cellID]=TileSpectra("maxADC",5,cellID,calib.GetTileCalib(cellID),event.GetROtype(),debug);
0741           hSpectraAlt[cellID].FillExtHGCROC(tempADC,toa,-1,nSampTOA,-1);
0742           hSpectraAlt[cellID].FillMaxVsTime(tempADC, aTile->GetCorrectedTOA(), nOffEmpty, nMaxADC);
0743         }
0744       }
0745     }
0746   }
0747 
0748   //==================================================================================
0749   // Setup general plotting infos
0750   //==================================================================================    
0751   
0752   TString outputDirPlots = GetPlotOutputDir();
0753   Double_t textSizeRel = 0.035;
0754 
0755   gSystem->Exec("mkdir -p "+outputDirPlots);
0756   
0757   StyleSettingsBasics("pdf");
0758   SetPlotStyle();  
0759   
0760   //==================================================================================
0761   // Create Summary plots
0762   //==================================================================================    
0763   TCanvas* canvas2DSigQA = new TCanvas("canvas2DSigQA","",0,0,1450,1300);  // gives the page size
0764   DefaultCanvasSettings( canvas2DSigQA, 0.08, 0.13, 0.045, 0.07);
0765 
0766   PlotSimple2D( canvas2DSigQA, hTOACorrVsCellID, -10000, setup->GetMaxCellID()+1, textSizeRel, Form("%s/TOACorrvsCellID.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
0767   PlotSimple2D( canvas2DSigQA, hTOACorrNsVsCellID, -10000, setup->GetMaxCellID()+1, textSizeRel, Form("%s/TOACorrNsvsCellID.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
0768   
0769   TCanvas* canvas1DSimple = new TCanvas("canvas1DSimple","",0,0,1450,1300);  // gives the page size
0770   DefaultCanvasSettings( canvas1DSimple, 0.08, 0.03, 0.03, 0.07);
0771 
0772   PlotSimple1D(canvas1DSimple, hSampleTOA, -10000, (double)it->second.samples, textSizeRel, Form("%s/NSampleToA.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1);
0773   PlotSimple1D(canvas1DSimple, hSampleMaxADC, -10000, (double)it->second.samples, textSizeRel, Form("%s/NSampleMaxADC.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1);
0774   PlotSimple1D(canvas1DSimple, hSampleAboveTh, -10000, (double)it->second.samples, textSizeRel, Form("%s/NSampleAboveTh.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1);
0775 
0776   // print general calib info
0777   calib.PrintGlobalInfo();  
0778 
0779   //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
0780   // Single layer plotting
0781   //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
0782   DetConf::Type detConf = setup->GetDetectorConfig();
0783   MultiCanvas panelPlot2D(detConf, "TimeWalk2D");
0784   bool init2D = panelPlot2D.Initialize(2);
0785     
0786   panelPlot2D.PlotCorr2DLayer(hSpectra, 1, -50, (it->second.samples)*25, -80, 800,
0787                               Form("%s/Waveform",outputDirPlots.Data()), plotSuffix.Data(), it->second, &calib );
0788   panelPlot2D.PlotCorr2DLayer(hSpectraAlt, 1, -50, (it->second.samples)*25, -80, 800,
0789                               Form("%s/WaveformMaxAdcOnly",outputDirPlots.Data()), plotSuffix.Data(), it->second, &calib );
0790   panelPlot2D.PlotCorr2DLayer(hSpectraAlt, 2, -6*1024, -2*1024, 0, 1000,
0791                               Form("%s/TOA_ADC_MaxAdcOnly",outputDirPlots.Data()), plotSuffix.Data(), it->second, &calib );
0792   
0793   //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
0794   // saving all outputs
0795   //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
0796   RootOutputHist->cd();
0797   if (typeRO == ReadOut::Type::Hgcroc){
0798     hSampleTOA->Write();
0799     hSampleMaxADC->Write();
0800     hSampleAboveTh->Write();
0801     hTOACorrVsCellID->Write();
0802     hTOACorrNsVsCellID->Write();
0803   }
0804 
0805   RootOutputHist->cd("IndividualCells");
0806   for(ithSpectra=hSpectra.begin(); ithSpectra!=hSpectra.end(); ++ithSpectra){
0807     ithSpectra->second.WriteExt(true);
0808   }
0809   RootOutputHist->cd("IndividualCellsAlt");
0810   for(ithSpectraAlt=hSpectraAlt.begin(); ithSpectraAlt!=hSpectraAlt.end(); ++ithSpectraAlt){
0811     ithSpectraAlt->second.WriteExt(true);
0812   }
0813 
0814   RootInput->Close();
0815   return true;
0816 }
0817 
0818 
0819 // ****************************************************************************
0820 // Investigate Cross Talk
0821 // ****************************************************************************
0822 bool HGCROC_Waveform_Analysis::InvestigateCrossTalk(void){
0823   std::cout<<"Analyse Waveform"<<std::endl;
0824   int evts=TdataIn->GetEntries();
0825 
0826   std::map<int,RunInfo> ri=readRunInfosFromFile(RunListInputName.Data(),debug,0);
0827   
0828   std::map<int,TileSpectra> hSpectraXtalk;
0829   std::map<int, TileSpectra>::iterator ithSpectraXtalk;
0830   std::map<int,TileSpectra> hSampleXtalk;
0831   std::map<int, TileSpectra>::iterator ithSampleXtalk;
0832   std::map<int,TileSpectra> hSpectraSat;
0833   std::map<int, TileSpectra>::iterator ithSpectraSat;
0834   std::map<int,TileSpectra> hSampleSat;
0835   std::map<int, TileSpectra>::iterator ithSampleSat;
0836   // set calib entry pointer
0837   TcalibIn->GetEntry(0);
0838 
0839   // obtain runnumber & RO-type
0840   TdataIn->GetEntry(0);
0841   Int_t runNr          = event.GetRunNumber();
0842   ReadOut::Type typeRO = event.GetROtype();
0843   std::map<int,RunInfo>::iterator it=ri.find(runNr);
0844 
0845   bool isSingleCh     = false;
0846   if (GetFixedROChannel() != -1) isSingleCh = true;
0847   
0848   int maxChannelPerLayer             = (setup->GetNMaxColumn()+1)*(setup->GetNMaxRow()+1)*(setup->GetNMaxModule()+1);
0849   int maxChannelPerLayerSingleMod    = (setup->GetNMaxColumn()+1)*(setup->GetNMaxRow()+1);
0850 
0851   
0852   //==================================================================================
0853   // create additional output hist
0854   //==================================================================================
0855   CreateOutputRootHistFile();
0856 
0857   TH1D* hNEvts                    = new TH1D("hNEvts","Events counter; ", 3, 0, 3);
0858   hNEvts->SetDirectory(0);
0859 
0860   TH2D* hSatADCVsLayer   = new TH2D( "hSatADCVsLayer","Saturated ADC; layer; brd channel; saturated/event ",
0861                                       setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
0862   hSatADCVsLayer->SetDirectory(0);
0863   hSatADCVsLayer->Sumw2();
0864   TH2D* hNegCellVsLayer   = new TH2D( "hNegCellVsLayer","Negative ADC; layer; brd channel; saturated/event ",
0865                                       setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
0866   hNegCellVsLayer->SetDirectory(0);
0867   hNegCellVsLayer->Sumw2();
0868 
0869   TH2D* h2DNSatCellsVsAsicHalfs = new TH2D( "h2DNSatCellsVsAsicHalfs","# Saturated per halfs; halfs; nCells sat. ",
0870                                             (setup->GetNMaxROUnit()+1)*2, 0, (setup->GetNMaxROUnit()+1)*2,
0871                                             36, -0.5, 36-0.5);
0872   h2DNSatCellsVsAsicHalfs->SetDirectory(0);
0873   h2DNSatCellsVsAsicHalfs->Sumw2();
0874   TH2D* h2DNNegCellsVsAsicHalfs = new TH2D( "h2DNNegCellsVsAsicHalfs","# Negative ADC per halfs; halfs; nCells w/ neg ADC ",
0875                                             (setup->GetNMaxROUnit()+1)*2, 0, (setup->GetNMaxROUnit()+1)*2,
0876                                             36, -0.5, 36-0.5);
0877   h2DNNegCellsVsAsicHalfs->SetDirectory(0);
0878   h2DNNegCellsVsAsicHalfs->Sumw2();
0879   TH2D* h2DNNegCellsVsAsicHalfsWT = new TH2D( "h2DNNegCellsVsAsicHalfsWT","# Negative ADC per halfs w/ ToA; halfs; nCells w/ neg ADC, w/ ToA ",
0880                                             (setup->GetNMaxROUnit()+1)*2, 0, (setup->GetNMaxROUnit()+1)*2,
0881                                             36, -0.5, 36-0.5);
0882   h2DNNegCellsVsAsicHalfsWT->SetDirectory(0);
0883   h2DNNegCellsVsAsicHalfsWT->Sumw2();
0884   
0885   TH2D* h2DNSatCellsVsNegCells = new TH2D( "h2DNSatCellsVsNegCells","# Saturated vs #Neg ADC; #cells sat; #cells neg ",
0886                                             100, -0.5, 100-0.5,
0887                                             setup->GetNActiveCells(), -0.5, setup->GetNActiveCells()-0.5);
0888   h2DNSatCellsVsNegCells->SetDirectory(0);
0889   h2DNSatCellsVsNegCells->Sumw2();
0890   TH2D* h2DNNegCellsVsNegCellsWToA = new TH2D( "h2DNNegCellsVsNegCellsWToA","# Neg ADC vs #Neg ADC w/ToA; #cells neg; #cells neg w/ ToA",
0891                                             setup->GetNActiveCells(), -0.5, setup->GetNActiveCells()-0.5,
0892                                             setup->GetNActiveCells(), -0.5, setup->GetNActiveCells()-0.5);
0893   h2DNNegCellsVsNegCellsWToA->SetDirectory(0);
0894   h2DNNegCellsVsNegCellsWToA->Sumw2();
0895   
0896   TH2D* h2DNSatCellsVsNegCellsWToA = new TH2D( "h2DNSatCellsVsNegCellsWToA","#Sat ADC vs #Neg ADC w/ToA; #cells sat; #cells neg w/ ToA ",
0897                                             100, -0.5, 100-0.5,
0898                                             setup->GetNActiveCells(), -0.5, setup->GetNActiveCells()-0.5);
0899   h2DNSatCellsVsNegCellsWToA->SetDirectory(0);
0900   h2DNSatCellsVsNegCellsWToA->Sumw2();
0901   
0902   TH2D* h2DAsicSatCellVsNegCell[setup->GetNMaxROUnit()+1][2];
0903   TH2D* h2DAsicSatCellVsNegCellWToA[setup->GetNMaxROUnit()+1][2];
0904   TH2D* h2DAsicNegCellVsNegCellWToA[setup->GetNMaxROUnit()+1][2];
0905   
0906   for (Int_t ro = 0; ro < setup->GetNMaxROUnit()+1; ro++){
0907     for (int h = 0; h< 2; h++){
0908       h2DAsicSatCellVsNegCell[ro][h]    = new TH2D( Form("h2DNSatCellsVsNegCells_%d_%d", ro, h),
0909                                                     "# Saturated vs #Neg ADC; #cells sat; #cells neg ",
0910                                                     32, -0.5, 32-0.5,
0911                                                     32, -0.5, 32-0.5);
0912       h2DAsicSatCellVsNegCell[ro][h]->SetDirectory(0);
0913       h2DAsicSatCellVsNegCell[ro][h]->Sumw2();
0914 
0915       h2DAsicSatCellVsNegCellWToA[ro][h]    = new TH2D( Form("h2DNSatCellsVsNegCellsWToA_%d_%d", ro, h),
0916                                                     "# Saturated vs #Neg ADC w/ TOA; #cells sat; #cells neg w/ ToA ",
0917                                                     32, -0.5, 32-0.5,
0918                                                     32, -0.5, 32-0.5);
0919       h2DAsicSatCellVsNegCellWToA[ro][h]->SetDirectory(0);
0920       h2DAsicSatCellVsNegCellWToA[ro][h]->Sumw2();
0921       h2DAsicNegCellVsNegCellWToA[ro][h]    = new TH2D( Form("h2DAsicNegCellVsNegCellWToA_%d_%d", ro, h),
0922                                                     "# Neg ADC vs #Neg ADC w/ TOA; #cells neg; #cells neg w/ ToA ",
0923                                                     32, -0.5, 32-0.5,
0924                                                     32, -0.5, 32-0.5);
0925       h2DAsicNegCellVsNegCellWToA[ro][h]->SetDirectory(0);
0926       h2DAsicNegCellVsNegCellWToA[ro][h]->Sumw2();
0927     }
0928   }
0929   
0930   RootOutputHist->mkdir("IndividualCellsXtalk");
0931   RootOutputHist->mkdir("IndividualCellsSat");
0932   RootOutputHist->cd("IndividualCellsXtalk");
0933   
0934   ROOT::Math::MinimizerOptions::SetDefaultMinimizer("Minuit2", "Migrad");  
0935   if (maxEvents == -1){
0936     maxEvents = evts;
0937   } else {
0938     std::cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl;
0939     std::cout << "ATTENTION: YOU ARE RESETTING THE MAXIMUM NUMBER OF EVENTS TO BE PROCESSED TO: " << maxEvents << ". THIS SHOULD ONLY BE USED FOR TESTING!" << std::endl;
0940     std::cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl;
0941   }
0942   
0943   //==================================================================================
0944   // process events from tree
0945   //==================================================================================
0946   for(int i=0; i<evts && i < maxEvents ; i++){
0947     if (i%5000 == 0&& i > 0 && debug > 0) std::cout << "Reading " <<  i << "/" << evts << " events"<< std::endl;
0948     if (debug > 2 ){
0949       std::cout << "************************************* NEW EVENT " << i << "  *********************************" << std::endl;
0950     }
0951     TdataIn->GetEntry(i);
0952     //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
0953     // process tiles in event to eval triggers
0954     //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++    
0955     int nSatTilesS  = 0;
0956     int nNegTilesS  = 0;
0957     int nSatTiles[(setup->GetNMaxROUnit()+1)*2];
0958     int nNegTiles[(setup->GetNMaxROUnit()+1)*2];
0959     int nToA          = 0;
0960     int nNegTilesSwT  = 0;
0961     int nNegTileswT[(setup->GetNMaxROUnit()+1)*2];
0962     for (int l = 0; l < (setup->GetNMaxROUnit()+1)*2; l++ ){
0963       nSatTiles[l]    = 0;
0964       nNegTiles[l]    = 0;
0965       nNegTileswT[l]  = 0;
0966     }
0967     
0968     for(int j=0; j<event.GetNTiles(); j++){
0969       Hgcroc* aTile=(Hgcroc*)event.GetTile(j);
0970       int cellID    = aTile->GetCellID();
0971       int layer     = setup->GetLayer(cellID);
0972       int chInLayer = setup->GetChannelInLayerFull(cellID);
0973       int asic      = setup->GetROunit(cellID);
0974       int roCh      = setup->GetROchannel(cellID);
0975       
0976       // Check if single channel eval is enabled
0977       if (isSingleCh){
0978         // only continue if we are looking at the channel under investigation
0979         if (GetFixedROChannel() != roCh)
0980           continue;
0981         else 
0982           if (debug> 3)std::cout << "triggered: " << (setup->DecodeCellID(cellID)).Data() << std::endl;
0983       }
0984       
0985       int aHalf     = int(roCh/36);
0986       int linHalfs  = asic*2+aHalf;
0987       double toa      = aTile->GetRawTOA();
0988       bool hasTOA     = false;
0989       // check if we have a toa
0990       if (toa > 0){
0991         hasTOA        = true;
0992         nToA++;
0993       }
0994       
0995       bool isSaturated = aTile->IsSaturatedADC();
0996       if (isSaturated){
0997         nSatTilesS++;
0998         nSatTiles[linHalfs]++;
0999         hSatADCVsLayer->Fill(layer,chInLayer);
1000       }
1001       int nSampBPed    = aTile->IsBelowPed(calib.GetPedestalSigH(cellID)*5.);
1002       bool isNegXtalk  = false;
1003       if (nSampBPed != -1){
1004         isNegXtalk     = true;
1005         nNegTilesS++;
1006         nNegTiles[linHalfs]++;
1007         hNegCellVsLayer->Fill(layer,chInLayer);
1008         if (hasTOA){
1009          nNegTilesSwT++;
1010          nNegTileswT[linHalfs]++;
1011         }
1012       }
1013     }
1014     hNEvts->Fill(0);
1015     if (nSatTilesS > 0){
1016       hNEvts->Fill(1);
1017       for (int l = 0; l < (setup->GetNMaxROUnit()+1)*2; l++ ){
1018         h2DNSatCellsVsAsicHalfs->Fill(l,nSatTiles[l]);
1019       }
1020     }
1021     if (nNegTilesS > 0){
1022       hNEvts->Fill(2);
1023       for (int l = 0; l < (setup->GetNMaxROUnit()+1)*2; l++ ){
1024         h2DNNegCellsVsAsicHalfs->Fill(l,nNegTiles[l]);
1025         h2DNNegCellsVsAsicHalfsWT->Fill(l,nNegTileswT[l]);
1026       }
1027     }
1028     
1029     h2DNSatCellsVsNegCells->Fill(nSatTilesS,nNegTilesS);
1030     h2DNSatCellsVsNegCellsWToA->Fill(nSatTilesS,nNegTilesSwT);
1031     h2DNNegCellsVsNegCellsWToA->Fill(nNegTilesS,nNegTilesSwT);
1032     
1033     for (int l = 0; l < (setup->GetNMaxROUnit()+1)*2; l++ ){
1034       int asic = int(l/2);
1035       int half = int(l%2);
1036       h2DAsicSatCellVsNegCell[asic][half]->Fill(nSatTiles[l],nNegTiles[l]);
1037       h2DAsicSatCellVsNegCellWToA[asic][half]->Fill(nSatTiles[l],nNegTileswT[l]);
1038       h2DAsicNegCellVsNegCellWToA[asic][half]->Fill(nNegTiles[l],nNegTileswT[l]);
1039     }
1040     
1041     if (nSatTilesS == 0 &&  nNegTilesS == 0)
1042       continue;
1043     else {
1044       if (debug > 3){
1045         if (nNegTilesS> 0){
1046           std::cout << "neg tile trigger: \t";
1047           for (int l = 0; l < (setup->GetNMaxROUnit()+1)*2; l++){
1048             std::cout << nNegTiles[l] << "\t";
1049           }
1050           std::cout << "\t tot\t: " << nNegTilesS <<  "\t" << "w/ toa\t" << nNegTilesSwT<< std::endl;
1051         }
1052         if (nSatTilesS >0){
1053           std::cout << "sat tile trigger: \t";
1054           for (int l = 0; l < (setup->GetNMaxROUnit()+1)*2; l++){
1055             std::cout << nSatTiles[l] << "\t";
1056           }
1057           std::cout << "\t tot\t: " << nSatTilesS << std::endl;
1058         }
1059       }
1060     }
1061     //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1062     // process tiles in event to fill selected events
1063     //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++    
1064     for(int j=0; j<event.GetNTiles(); j++){
1065       //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++  
1066       // histo filling for HGCROC specific things & resetting of ped
1067       //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++        
1068       Hgcroc* aTile=(Hgcroc*)event.GetTile(j);
1069       int cellID      = aTile->GetCellID();
1070       int roCh        = setup->GetROchannel(cellID);
1071       int aHalf       = int(roCh/36);
1072       int asic        = setup->GetROunit(cellID);
1073       int linHalfs    = asic*2+aHalf;
1074 
1075       ithSpectraSat   = hSpectraSat.find(cellID);
1076       ithSpectraXtalk = hSpectraXtalk.find(cellID);
1077       ithSampleSat    = hSampleSat.find(cellID);
1078       ithSampleXtalk  = hSampleXtalk.find(cellID);
1079       
1080       // get pedestal values from calib object
1081       double ped    = aTile->GetPedestal();
1082       
1083       // obtain integrated tile quantities for histo filling
1084       double toaRaw   = aTile->GetRawTOA();
1085       bool hasTOA     = false;
1086       Int_t nOffEmpty = 2;
1087       
1088       // check if we have a toa
1089       if (toaRaw > 0)
1090         hasTOA        = true;
1091 
1092       int trigAsicH = 0;
1093       if (GetFixedROChannel() > 31)
1094         trigAsicH   = 1;
1095       
1096       // Detailed debug info with print of waveform
1097       if (debug > 3){
1098         aTile->PrintWaveFormDebugInfo(calib.GetPedestalMeanH(cellID), calib.GetPedestalMeanL(cellID), calib.GetPedestalSigL(cellID));
1099       }
1100     
1101       //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1102       // fill as function of sample, no TOA requirement
1103       //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1104       if (!isSingleCh){
1105         if (nSatTilesS > 0){
1106           if(ithSampleSat!=hSampleSat.end()){
1107             ithSampleSat->second.FillWaveform(aTile->GetADCWaveform(), ped);
1108           } else {
1109             RootOutputHist->cd("IndividualCellsSat");
1110             hSampleSat[cellID]=TileSpectra("SatSample",7,cellID,calib.GetTileCalib(cellID),event.GetROtype(),debug);
1111             hSampleSat[cellID].FillWaveform(aTile->GetADCWaveform(), ped);
1112           }
1113         }
1114         if (nNegTilesS > 0 ){
1115           if(ithSampleXtalk!=hSampleXtalk.end()){
1116             ithSampleXtalk->second.FillWaveform(aTile->GetADCWaveform(), ped);
1117           } else {
1118             RootOutputHist->cd("IndividualCellsXtalk");
1119             hSampleXtalk[cellID]=TileSpectra("XtalkSample",7,cellID,calib.GetTileCalib(cellID),event.GetROtype(),debug);
1120             hSampleXtalk[cellID].FillWaveform(aTile->GetADCWaveform(), ped);
1121           }
1122         }
1123       } else {
1124         if (nSatTiles[asic*2 + trigAsicH] > 0 ){
1125           if(ithSampleSat!=hSampleSat.end()){
1126             ithSampleSat->second.FillWaveform(aTile->GetADCWaveform(), ped);
1127           } else {
1128             RootOutputHist->cd("IndividualCellsSat");
1129             hSampleSat[cellID]=TileSpectra("SatSample",7,cellID,calib.GetTileCalib(cellID),event.GetROtype(),debug);
1130             hSampleSat[cellID].FillWaveform(aTile->GetADCWaveform(), ped);
1131           }
1132         }
1133         
1134         if (nNegTiles[asic*2 + trigAsicH] > 0 ){
1135           if(ithSampleXtalk!=hSampleXtalk.end()){
1136             ithSampleXtalk->second.FillWaveform(aTile->GetADCWaveform(), ped);
1137           } else {
1138             RootOutputHist->cd("IndividualCellsXtalk");
1139             hSampleXtalk[cellID]=TileSpectra("XtalkSample",7,cellID,calib.GetTileCalib(cellID),event.GetROtype(),debug);
1140             hSampleXtalk[cellID].FillWaveform(aTile->GetADCWaveform(), ped);
1141           }
1142         }         
1143       }
1144       
1145     
1146       // check whether tile has TOA, don't continue if it doesn't
1147       if (!hasTOA) continue;
1148     
1149       //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1150       // fill as function of time, needs TOA
1151       //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1152       if (!isSingleCh){
1153         if (nSatTilesS > 0){
1154           if(ithSpectraSat!=hSpectraSat.end()){
1155             ithSpectraSat->second.FillWaveformVsTime(aTile->GetADCWaveform(), aTile->GetCorrectedTOA(), ped, nOffEmpty);
1156           } else {
1157             RootOutputHist->cd("IndividualCellsSat");
1158             hSpectraSat[cellID]=TileSpectra("Saturated",6,cellID,calib.GetTileCalib(cellID),event.GetROtype(),debug);
1159             hSpectraSat[cellID].FillWaveformVsTime(aTile->GetADCWaveform(), aTile->GetCorrectedTOA(), ped, nOffEmpty);
1160           }
1161           
1162         }
1163         if (nNegTilesS > 0 ){
1164           if(ithSpectraXtalk!=hSpectraXtalk.end()){
1165             ithSpectraXtalk->second.FillWaveformVsTime(aTile->GetADCWaveform(), aTile->GetCorrectedTOA(), ped, nOffEmpty);
1166           } else {
1167             RootOutputHist->cd("IndividualCellsXtalk");
1168             hSpectraXtalk[cellID]=TileSpectra("Xtalk",6,cellID,calib.GetTileCalib(cellID),event.GetROtype(),debug);
1169             hSpectraXtalk[cellID].FillWaveformVsTime(aTile->GetADCWaveform(), aTile->GetCorrectedTOA(), ped, nOffEmpty);
1170           }
1171         }
1172       } else {
1173         if (nSatTiles[asic*2 + trigAsicH] > 0 ){
1174           if(ithSpectraSat!=hSpectraSat.end()){
1175             ithSpectraSat->second.FillWaveformVsTime(aTile->GetADCWaveform(), aTile->GetCorrectedTOA(), ped, nOffEmpty);
1176           } else {
1177             RootOutputHist->cd("IndividualCellsSat");
1178             hSpectraSat[cellID]=TileSpectra("Saturated",6,cellID,calib.GetTileCalib(cellID),event.GetROtype(),debug);
1179             hSpectraSat[cellID].FillWaveformVsTime(aTile->GetADCWaveform(), aTile->GetCorrectedTOA(), ped, nOffEmpty);
1180           }
1181         }
1182         
1183         if (nNegTiles[asic*2 + trigAsicH] > 0 ){
1184           if(ithSpectraXtalk!=hSpectraXtalk.end()){
1185             ithSpectraXtalk->second.FillWaveformVsTime(aTile->GetADCWaveform(), aTile->GetCorrectedTOA(), ped, nOffEmpty);
1186           } else {
1187             RootOutputHist->cd("IndividualCellsXtalk");
1188             hSpectraXtalk[cellID]=TileSpectra("Xtalk",6,cellID,calib.GetTileCalib(cellID),event.GetROtype(),debug);
1189             hSpectraXtalk[cellID].FillWaveformVsTime(aTile->GetADCWaveform(), aTile->GetCorrectedTOA(), ped, nOffEmpty);
1190           }
1191         } 
1192       }
1193     }
1194   }
1195   
1196   RootInput->Close();
1197 
1198   //==================================================================================
1199   // write outputs to root fiif (de)le
1200   //==================================================================================
1201   RootOutputHist->cd();
1202     h2DNNegCellsVsAsicHalfs->Write();
1203     h2DNNegCellsVsAsicHalfsWT->Write();
1204     h2DNSatCellsVsAsicHalfs->Write();
1205     h2DNSatCellsVsNegCells->Write();
1206     h2DNSatCellsVsNegCellsWToA->Write();
1207     h2DNNegCellsVsNegCellsWToA->Write();
1208     hSatADCVsLayer->Write();
1209     hNegCellVsLayer->Write();
1210     
1211     for (Int_t ro = 0; ro < setup->GetNMaxROUnit()+1; ro++){
1212       for (int h = 0; h< 2; h++){
1213         h2DAsicSatCellVsNegCell[ro][h]->Write();
1214         h2DAsicSatCellVsNegCellWToA[ro][h]->Write();
1215         h2DAsicNegCellVsNegCellWToA[ro][h]->Write();
1216       }
1217     }
1218     hNEvts->Write();
1219 
1220 
1221   RootOutputHist->cd("IndividualCellsXtalk");
1222     for(ithSpectraXtalk=hSpectraXtalk.begin(); ithSpectraXtalk!=hSpectraXtalk.end(); ++ithSpectraXtalk){
1223       ithSpectraXtalk->second.WriteExt(true);
1224     }
1225     for(ithSampleXtalk=hSampleXtalk.begin(); ithSampleXtalk!=hSampleXtalk.end(); ++ithSampleXtalk){
1226       ithSampleXtalk->second.WriteExt(true);
1227     }
1228   RootOutputHist->cd("IndividualCellsSat");
1229     for(ithSpectraSat=hSpectraSat.begin(); ithSpectraSat!=hSpectraSat.end(); ++ithSpectraSat){
1230       ithSpectraSat->second.WriteExt(true);
1231     }
1232     for(ithSampleSat=hSampleSat.begin(); ithSampleSat!=hSampleSat.end(); ++ithSampleSat){
1233       ithSampleSat->second.WriteExt(true);
1234     }
1235   RootOutputHist->Close();
1236 
1237   // create histos if they don't exist for plotting reasons
1238   for (int id = 0; id < setup->GetMaxCellID()+1; id++){
1239     // skip if not defined in setup
1240     if (!setup->ContainedInSetup(id)) continue; 
1241     ithSpectraSat   = hSpectraSat.find(id);
1242     ithSpectraXtalk = hSpectraXtalk.find(id);
1243     ithSampleSat    = hSampleSat.find(id);
1244     ithSampleXtalk  = hSampleXtalk.find(id);
1245     
1246     if(ithSpectraXtalk==hSpectraXtalk.end()){
1247       hSpectraXtalk[id]=TileSpectra("Xtalk",6,id,calib.GetTileCalib(id),typeRO,debug);
1248     }
1249     if(ithSpectraSat==hSpectraSat.end()){
1250       hSpectraSat[id]=TileSpectra("Saturated",6,id,calib.GetTileCalib(id),typeRO,debug);
1251     }
1252     if(ithSampleXtalk==hSampleXtalk.end()){
1253       hSampleXtalk[id]=TileSpectra("XtalkSample",7,id,calib.GetTileCalib(id),typeRO,debug);
1254     }
1255     if(ithSampleSat==hSampleSat.end()){
1256       hSampleSat[id]=TileSpectra("SatSample",7,id,calib.GetTileCalib(id),typeRO,debug);
1257     }
1258   }
1259   
1260   //==================================================================================
1261   // Setup general plotting infos
1262   //==================================================================================    
1263   
1264   TString outputDirPlots = GetPlotOutputDir();
1265   Double_t textSizeRel = 0.035;
1266 
1267   if (isSingleCh){
1268     outputDirPlots = Form("%s/Channel_%d", outputDirPlots.Data(), GetFixedROChannel());
1269   }
1270   gSystem->Exec("mkdir -p "+outputDirPlots);
1271 
1272   StyleSettingsBasics("pdf");
1273   SetPlotStyle();  
1274   
1275   //==================================================================================
1276   // Create Summary plots
1277   //==================================================================================    
1278   TCanvas* canvas2DSigQA = new TCanvas("canvas2DSigQA","",0,0,1450,1300);  // gives the page size
1279   DefaultCanvasSettings( canvas2DSigQA, 0.08, 0.13, 0.045, 0.07);
1280 
1281   canvas2DSigQA->SetLogz();
1282   
1283   h2DNNegCellsVsAsicHalfs->Scale(1./maxEvents);
1284   PlotSimple2D( canvas2DSigQA, h2DNNegCellsVsAsicHalfs, -10000, -10000, textSizeRel, Form("%s/NegCellsvsAsicHalfs.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);  
1285   h2DNNegCellsVsAsicHalfsWT->Scale(1./maxEvents);
1286   PlotSimple2D( canvas2DSigQA, h2DNNegCellsVsAsicHalfsWT, -10000, -10000, textSizeRel, Form("%s/NegCellsvsAsicHalfsWithToA.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
1287   h2DNSatCellsVsAsicHalfs->Scale(1./maxEvents);
1288   PlotSimple2D( canvas2DSigQA, h2DNSatCellsVsAsicHalfs, 20, -10000, textSizeRel, Form("%s/SatCellsvsAsicHalfs.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
1289 
1290   h2DNSatCellsVsNegCells->Scale(1./maxEvents);
1291   PlotSimple2D( canvas2DSigQA, h2DNSatCellsVsNegCells, 200, 20, textSizeRel, Form("%s/SatCellsvsNegCells.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
1292   h2DNSatCellsVsNegCellsWToA->Scale(1./maxEvents);
1293   PlotSimple2D( canvas2DSigQA, h2DNSatCellsVsNegCellsWToA, 100, 20, textSizeRel, Form("%s/SatCellsvsNegCellsWToA.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
1294   h2DNNegCellsVsNegCellsWToA->Scale(1./maxEvents);
1295   PlotSimple2D( canvas2DSigQA, h2DNNegCellsVsNegCellsWToA, 20, 200, textSizeRel, Form("%s/NegCellsvsNegCellsWToA.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
1296  
1297   hSatADCVsLayer->Scale(1./maxEvents);
1298   PlotSimple2D( canvas2DSigQA, hSatADCVsLayer, -10000, -10000, textSizeRel, Form("%s/SatCellsvsLayer.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
1299   hNegCellVsLayer->Scale(1./maxEvents);
1300   PlotSimple2D( canvas2DSigQA, hNegCellVsLayer, -10000, -10000, textSizeRel, Form("%s/NegCellsvsLayer.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);  
1301   
1302   for (Int_t ro = 0; ro < setup->GetNMaxROUnit()+1; ro++){
1303     for (int h = 0; h< 2; h++){
1304       h2DAsicSatCellVsNegCell[ro][h]->Scale(1./maxEvents);
1305       PlotSimple2D( canvas2DSigQA, h2DAsicSatCellVsNegCell[ro][h], -10000, 10,  textSizeRel,
1306                         Form("%s/SatCellVsNegCell_Asic_%d_Half_%d.%s", outputDirPlots.Data(), ro, h, plotSuffix.Data()), 
1307                         it->second, 1, kFALSE, "colz",true, Form("Asic %d, Half %d",ro,h));
1308       h2DAsicSatCellVsNegCellWToA[ro][h]->Scale(1./maxEvents);
1309       PlotSimple2D( canvas2DSigQA, h2DAsicSatCellVsNegCellWToA[ro][h], -10000, 10, textSizeRel,
1310                         Form("%s/SatCellsvsNegCellsWToA_Asic_%d_Half_%d.%s", outputDirPlots.Data(), ro, h, plotSuffix.Data()), 
1311                         it->second, 1, kFALSE, "colz",true, Form("Asic %d, Half %d",ro,h));
1312       h2DAsicNegCellVsNegCellWToA[ro][h]->Scale(1./maxEvents);
1313       PlotSimple2D( canvas2DSigQA, h2DAsicNegCellVsNegCellWToA[ro][h], -10000, -10000, textSizeRel,
1314                         Form("%s/NegCellsvsNegCellsWToA_Asic_%d_Half_%d.%s", outputDirPlots.Data(), ro, h, plotSuffix.Data()), 
1315                         it->second, 1, kFALSE, "colz",true, Form("Asic %d, Half %d",ro,h));
1316     }
1317   }
1318   
1319   TCanvas* canvas1DSimple = new TCanvas("canvas1DSimple","",0,0,1450,1300);  // gives the page size
1320   DefaultCanvasSettings( canvas1DSimple, 0.08, 0.03, 0.03, 0.07);
1321   
1322   PlotSimple1D(canvas1DSimple, hNEvts, -10000, -10000, textSizeRel, Form("%s/NEvts.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1);
1323 
1324   // print general calib info
1325   calib.PrintGlobalInfo();  
1326 
1327   //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++  
1328   // Single Layer Plotting
1329   //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++  
1330   DetConf::Type detConf = setup->GetDetectorConfig();
1331   MultiCanvas panelPlot2D(detConf, "Saturation");
1332   bool init2D = panelPlot2D.Initialize(2);
1333   
1334   if (ExtPlot > 1 ){
1335     std::cout << "plotting layers" << std::endl;
1336     panelPlot2D.PlotCorr2DLayer(hSpectraXtalk, 1, -50, (it->second.samples)*25, -80, 800,
1337                                 Form("%s/WaveformNeg",outputDirPlots.Data()), plotSuffix.Data(), it->second, &calib );
1338     panelPlot2D.PlotCorr2DLayer(hSpectraSat, 1, -50, (it->second.samples)*25, -80, 800,
1339                                 Form("%s/WaveformSat",outputDirPlots.Data()), plotSuffix.Data(), it->second, &calib );
1340   }
1341 
1342   //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++  
1343   // Single Asic Plotting
1344   //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++  
1345   std::cout << "plotting full asic" << std::endl;  
1346   MultiCanvas panelPlotAsic(DetConf::Type::Asic, "Saturation");
1347   bool initAsic = panelPlotAsic.Initialize(2);
1348 
1349   panelPlotAsic.PlotCorr2DLayer(hSpectraXtalk, 1, -50, (it->second.samples)*25, -80, 800,
1350                                 Form("%s/WaveformNeg",outputDirPlots.Data()), plotSuffix.Data(), it->second, &calib, 0, GetFixedROChannel() );
1351   panelPlotAsic.PlotCorr2DLayer(hSpectraSat, 1, -50, (it->second.samples)*25, -80, 800,
1352                                 Form("%s/WaveformSat",outputDirPlots.Data()), plotSuffix.Data(), it->second, &calib, 0, GetFixedROChannel() );
1353   panelPlotAsic.PlotCorr2DLayer(hSampleXtalk, 1, 0, (it->second.samples), -80, 800,
1354                                 Form("%s/WaveformSampleNeg",outputDirPlots.Data()), plotSuffix.Data(), it->second, &calib, 0, GetFixedROChannel() );
1355   panelPlotAsic.PlotCorr2DLayer(hSampleSat, 1, 0, (it->second.samples), -80, 800,
1356                                 Form("%s/WaveformSampleSat",outputDirPlots.Data()), plotSuffix.Data(), it->second, &calib, 0, GetFixedROChannel() );
1357   std::cout << "ending plotting full asic" << std::endl;    
1358   return true;
1359 }
1360 
1361 
1362 //***********************************************************************************************
1363 //*********************** Create output files ***************************************************
1364 //***********************************************************************************************
1365 bool HGCROC_Waveform_Analysis::CreateOutputRootFile(void){
1366   if(Overwrite){
1367     std::cout << "overwriting exisiting output file" << std::endl;
1368     RootOutput=new TFile(RootOutputName.Data(),"RECREATE");
1369   } else{
1370     std::cout << "creating output file" << std::endl;
1371     RootOutput = new TFile(RootOutputName.Data(),"CREATE");
1372   }
1373   if(RootOutput->IsZombie()){
1374     std::cout<<"Error opening '"<<RootOutput<<"'no reachable path? Exist without force mode to overwrite?..."<<std::endl;
1375     return false;
1376   }
1377   return true;
1378 }
1379 
1380 //***********************************************************************************************
1381 //*********************** Create output files ***************************************************
1382 //***********************************************************************************************
1383 bool HGCROC_Waveform_Analysis::CreateOutputRootHistFile(void){
1384   std::cout << "Additional Output with histos being created: " << RootOutputNameHist.Data() << std::endl;
1385   if(Overwrite){
1386     std::cout << "overwriting exisiting output file" << std::endl;
1387     RootOutputHist=new TFile(RootOutputNameHist.Data(),"RECREATE");
1388   } else{
1389     std::cout << "creating output file" << std::endl;
1390     RootOutputHist = new TFile(RootOutputNameHist.Data(),"CREATE");
1391   }
1392   if(RootOutputHist->IsZombie()){
1393     std::cout<<"Error opening '"<<RootOutputHist<<"'no reachable path? Exist without force mode to overwrite?..."<<std::endl;
1394     return false;
1395   }
1396   return true;
1397 }
[ 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!