Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-09-15 08:17:56 

0001 #include "Analyses.h"
0002 #include <vector>
0003 #include "TROOT.h"
0004 #include <unistd.h>
0005 #include "TF1.h"
0006 #include "TFitResult.h"
0007 #include "TFitResultPtr.h"
0008 #include "TH1D.h"
0009 #include "TH2D.h"
0010 #include "TProfile.h"
0011 #include "TChain.h"
0012 #include "CommonHelperFunctions.h"
0013 #include "TileSpectra.h"
0014 #include "PlottHelper.h"
0015 #include "TRandom3.h"
0016 #include "TMath.h"
0017 #include "Math/Minimizer.h"
0018 #include "Math/MinimizerOptions.h"
0019 
0020 #include "HGCROC_Convert.h"
0021 
0022 #include "waveform_fitting/waveform_fit_base.h"
0023 #include "waveform_fitting/crystal_ball_fit.h"
0024 #include "waveform_fitting/max_sample_fit.h"
0025 
0026 // ****************************************************************************
0027 // Checking and opening input and output files
0028 // ****************************************************************************
0029 bool Analyses::CheckAndOpenIO(void){
0030   int matchingbranch;
0031   if(!ASCIIinputName.IsNull()){
0032     std::cout << "Input to be converted into correct format :" <<  ASCIIinputName.Data() << std::endl;
0033     ASCIIinput.open(ASCIIinputName.Data(),std::ios::in);
0034     if(!ASCIIinput.is_open()){
0035       std::cout<<"Could not open input file: "<<std::endl;
0036       return false;
0037     }
0038   }
0039 
0040   //Need to check first input to get the setup...I do not think it is necessary
0041   std::cout <<"=============================================================" << std::endl;
0042   std::cout<<"Input name set to: "<<RootInputName.Data() <<std::endl;
0043   std::cout<<"Output name set to: "<<RootOutputName.Data() <<std::endl;
0044   if (!Convert) std::cout<<"Calib name set to: "<<RootCalibInputName.Data() <<std::endl;
0045   std::cout <<"=============================================================" << std::endl;
0046   if(!RootInputName.IsNull()){
0047     //File exist?
0048     RootInput=new TFile(RootInputName.Data(),"READ");
0049     if(RootInput->IsZombie()){
0050       std::cout<<"Error opening '"<<RootInputName<<"', does the file exist?"<<std::endl;
0051       return false;
0052     }
0053 
0054     //Retrieve info, start with setup
0055     TsetupIn = (TTree*) RootInput->Get("Setup");
0056     if(!TsetupIn){
0057       std::cout<<"Could not retrieve the Setup tree, leaving"<<std::endl;
0058       return false;
0059     }
0060     setup=Setup::GetInstance();
0061     std::cout<<"Setup add "<<setup<<std::endl;
0062     //matchingbranch=TsetupIn->SetBranchAddress("setup",&setup);
0063     matchingbranch=TsetupIn->SetBranchAddress("setup",&rswptr);
0064     if(matchingbranch<0){
0065       std::cout<<"Error retrieving Setup info from the tree"<<std::endl;
0066       return false;
0067     }
0068     std::cout<<"Entries "<<TsetupIn->GetEntries()<<std::endl;
0069     TsetupIn->GetEntry(0);
0070     setup->Initialize(*rswptr);
0071     std::cout<<"Reading "<<RootInput->GetName()<<std::endl;
0072     std::cout<<"Setup Info "<<setup->IsInit()<<"  and  "<<setup->GetCellID(0,0)<<std::endl;
0073     //std::cout<<"Setup add now "<<setup<<std::endl;
0074 
0075     //Retrieve info, existing data?
0076     TdataIn = (TTree*) RootInput->Get("Data");
0077     if(!TdataIn){
0078       std::cout<<"Could not retrieve the Data tree, leaving"<<std::endl;
0079       return false;
0080     }
0081     matchingbranch=TdataIn->SetBranchAddress("event",&eventptr);
0082     if(matchingbranch<0){
0083       std::cout<<"Error retrieving Event info from the tree"<<std::endl;
0084       return false;
0085     }
0086     
0087     //Do I really want this?
0088     TcalibIn = (TTree*) RootInput->Get("Calib");
0089     if(!TcalibIn){
0090       std::cout<<"Could not retrieve Calib tree, leaving"<<std::endl;
0091       //return false;
0092     }
0093     else {
0094       matchingbranch=TcalibIn->SetBranchAddress("calib",&calibptr);
0095       if(matchingbranch<0){
0096         std::cout<<"Error retrieving calibration info from the tree"<<std::endl;
0097         TcalibIn=nullptr;
0098       }
0099     }
0100     //End of do I really want this?
0101     
0102     //We want to retrieve also calibration if do not specify ApplyTransferCalib from external file
0103     //In other words, the pedestal was potentially already done and we have an existing calib object
0104     if((!ApplyTransferCalib && ExtractScaling) || (!ApplyTransferCalib && ExtractScalingImproved) || (!ApplyTransferCalib && ReextractNoise)){
0105       TcalibIn = (TTree*) RootInput->Get("Calib");
0106       if(!TcalibIn){
0107         std::cout<<"Could not retrieve Calib tree, leaving"<<std::endl;
0108         return false;
0109       }
0110       //calib=Calib::GetInstance();
0111       matchingbranch=TcalibIn->SetBranchAddress("calib",&calibptr);
0112       if(matchingbranch<0){
0113         std::cout<<"Error retrieving calibration info from the tree"<<std::endl;
0114         return false;
0115       }
0116     }
0117     //All good
0118   }
0119   else if(!Convert){
0120     std::cout<<"Explicit Input option mandatory except for convertion ASCII -> ROOT"<<std::endl;
0121     return false;
0122   }  
0123   
0124   //Setup Output files
0125   if(!RootOutputName.IsNull()){    
0126     if (!Convert){
0127       TString temp = RootOutputName;
0128       temp         = temp.ReplaceAll(".root","_Hists.root");
0129       SetRootOutputHists(temp);
0130       // std::cout << "creating additional histo file: " << RootOutputNameHist.Data() << " tree in : "<< RootOutputName.Data() << std::endl;
0131     }
0132     
0133     bool sCOF = CreateOutputRootFile();
0134     if (!sCOF) return false;
0135     
0136     TsetupOut = new TTree("Setup","Setup");
0137     setup=Setup::GetInstance();
0138     //TsetupOut->Branch("setup",&setup);
0139     TsetupOut->Branch("setup",&rsw);
0140     TdataOut = new TTree("Data","Data");
0141     TdataOut->Branch("event",&event);
0142     TcalibOut = new TTree("Calib","Calib");
0143     TcalibOut->Branch("calib",&calib);
0144   }
0145   else if (!SaveCalibOnly){
0146     std::cout<<"Output option mandatory except when converting"<<std::endl;
0147     return false;
0148   }
0149   
0150   if(!RootCalibInputName.IsNull()){
0151     RootCalibInput=new TFile(RootCalibInputName.Data(),"READ");
0152     if(RootCalibInput->IsZombie()){
0153       std::cout<<"Error opening '"<<RootCalibInputName<<"', does the file exist?"<<std::endl;
0154       return false;
0155     }
0156     TcalibIn = nullptr;
0157     TcalibIn = (TTree*) RootCalibInput->Get("Calib");
0158     if(!TcalibIn){
0159       std::cout<<"Could not retrieve Calib tree, leaving"<<std::endl;
0160       return false;
0161     } else {
0162       std::cout<<"Retrieved calib object from external input!"<<std::endl;
0163     }
0164     matchingbranch=TcalibIn->SetBranchAddress("calib",&calibptr);
0165     if(matchingbranch<0){
0166       std::cout<<"Error retrieving calibration info from the tree"<<std::endl;
0167       return false;
0168     }else {
0169       std::cout<<"Correctly set branch for external calib input."<<std::endl;
0170     }
0171     
0172   }
0173 
0174   if(!RootCalibOutputName.IsNull() && SaveCalibOnly){
0175     std::cout << "entered here" << std::endl;
0176     RootCalibOutput=new TFile(RootCalibOutputName.Data(),"RECREATE");
0177     if(RootCalibOutput->IsZombie()){
0178       std::cout<<"Error opening '"<<RootCalibOutputName<<"', does the file exist?"<<std::endl;
0179       return false;
0180     }
0181     
0182     if (RootOutputName.IsNull()){
0183       TsetupOut = new TTree("Setup","Setup");
0184       setup=Setup::GetInstance();
0185       //TsetupOut->Branch("setup",&setup);
0186       TsetupOut->Branch("setup",&rsw);
0187       TcalibOut = new TTree("Calib","Calib");
0188       TcalibOut->Branch("calib",&calib);
0189     }
0190   }
0191 
0192   if(!RootPedestalInputName.IsNull()){
0193     RootPedestalInput = new TFile(RootPedestalInputName.Data(),"READ");
0194     if(RootPedestalInput->IsZombie()){
0195       std::cout<<"Error opening '"<<RootPedestalInputName<<"', does the file exist?"<<std::endl;
0196       return false;
0197     }
0198     TcalibIn = (TTree*) RootPedestalInput->Get("Calib");
0199     if(!TcalibIn){
0200       std::cout<<"Could not retrieve Calib tree, leaving"<<std::endl;
0201       return false;
0202     } else {
0203       std::cout<<"Retrieved calib object from external input for pedestals!"<<std::endl;
0204     }
0205     matchingbranch=TcalibIn->SetBranchAddress("calib",&calibptr);
0206     if(matchingbranch<0){
0207       std::cout<<"Error retrieving calibration info from the tree"<<std::endl;
0208       return false;
0209     }else {
0210       std::cout<<"Correctly set branch for external calib for pedestal input."<<std::endl;
0211     }
0212     //std::cout<<"Did the address changed? "<<&calib<<std::endl;
0213   }
0214   if(!MapInputName.IsNull()){
0215     MapInput.open(MapInputName.Data(),std::ios::in);
0216     if(!MapInput.is_open()){
0217       std::cout<<"Could not open mapping file: "<<MapInputName<<std::endl;
0218       return false;
0219     }   
0220   }
0221   std::cout <<"=============================================================" << std::endl;
0222   std::cout <<" Basic setup complete" << std::endl;
0223   std::cout <<"=============================================================" << std::endl;
0224   return true;    
0225 }
0226 
0227 // ****************************************************************************
0228 // Primary process function to start all calibrations
0229 // ****************************************************************************
0230 bool Analyses::Process(void){
0231   bool status = true;
0232   ROOT::EnableImplicitMT();
0233   
0234   if(Convert){
0235     std::cout << "Converting !" << std::endl;
0236     if (HGCROC) {
0237       // Set waveform builder
0238       waveform_fit_base *waveform_builder = nullptr;
0239       // waveform_builder = new crystal_ball_fit();
0240 
0241       // // Set the parameters based off what I found in the stand alone analysis
0242       // // Alpha
0243       // waveform_builder->set_parameter(0, 1.1);  // Initial value
0244       // waveform_builder->set_parameter(10, 1);   // Lower bound
0245       // waveform_builder->set_parameter(20, 1.2); // Uppser bound
0246 
0247       // // n
0248       // waveform_builder->set_parameter(1, 0.3);
0249       // waveform_builder->set_parameter(11, 0.2);
0250       // waveform_builder->set_parameter(21, 0.6);
0251 
0252       // // x_bar
0253       // waveform_builder->set_parameter(2, 0.7);
0254       // waveform_builder->set_parameter(12, 0.5);
0255       // waveform_builder->set_parameter(22, 4.5);
0256 
0257       // // sigma
0258       // waveform_builder->set_parameter(3, 0.3);
0259       // waveform_builder->set_parameter(13, 0.25);
0260       // waveform_builder->set_parameter(23, 0.65);
0261 
0262       // // N
0263       // waveform_builder->set_parameter(4, 0);
0264       // waveform_builder->set_parameter(14, 0);
0265       // waveform_builder->set_parameter(24, 2000);
0266 
0267       // // Offset
0268       // waveform_builder->set_parameter(5, 100);  // This needs to become pedestals eventually... 
0269       // waveform_builder->set_parameter(15, 0);
0270       // waveform_builder->set_parameter(25, 160);
0271 
0272       waveform_builder = new max_sample_fit();
0273 
0274       std::cout << "Running HGCROC conversion" << std::endl;
0275       status=run_hgcroc_conversion(this, waveform_builder);
0276     } else {
0277       if (!(GetASCIIinputName().EndsWith(".root"))){
0278         status=ConvertASCII2Root();
0279       } else {
0280         std::cout << "WARNING: This option should only be used for the 2023 SPS test beam for which the CAEN raw data was lost!" << std::endl;
0281         status=ConvertOldRootFile2Root();
0282       }
0283     }
0284   }
0285   
0286   // extract pedestal from pure pedestal runs (internal triggers)
0287   if(ExtractPedestal){
0288     status=GetPedestal();
0289   }
0290   
0291   // copy existing calibration to other file
0292   if(ApplyTransferCalib){
0293     status=TransferCalib();
0294   }
0295   
0296   // extract mip calibration 
0297   if(ExtractScaling){
0298     status=GetScaling();
0299   }
0300   
0301   // extract noise sample from trigger flagged files
0302   if(ReextractNoise){
0303     status=GetNoiseSampleAndRefitPedestal();
0304   }
0305   
0306   // rerun mip calibration based on triggers
0307   if (ExtractScalingImproved){
0308     status=GetImprovedScaling();
0309   }
0310   
0311   // copy full calibration to different file and calculate energy
0312   if(ApplyCalibration){
0313     status=Calibrate();
0314   }
0315   
0316   // reduce file to only noise triggers
0317   if(SaveNoiseOnly){
0318     status=SaveNoiseTriggersOnly();
0319   }
0320   
0321   // reduce file to only mip triggers
0322   if(SaveMipsOnly){
0323     status=SaveMuonTriggersOnly();
0324   }
0325   
0326   // skim HGCROC data to discard pure noise events
0327   if(SkimHGCROC){
0328     status=SkimHGCROCData();
0329   }
0330   
0331   // reduce file to only mip triggers
0332   if(EvalLocalTriggers){
0333     status=RunEvalLocalTriggers();
0334   }
0335   
0336   // save calib to output only
0337   if (SaveCalibOnly){
0338     status = SaveCalibToOutputOnly();
0339   }
0340   return status;
0341 }
0342 
0343 
0344 // ****************************************************************************
0345 // convert original ASCII file from CAEN output to root format
0346 // ****************************************************************************
0347 bool Analyses::ConvertASCII2Root(void){
0348   //============================================
0349   //Init first
0350   //============================================
0351   // initialize setup
0352   if (MapInputName.CompareTo("")== 0) {
0353       std::cerr << "ERROR: No mapping file has been provided, please make sure you do so! Aborting!" << std::endl;
0354       return false;
0355   }
0356   std::cout << "creating mapping " << std::endl;
0357   setup->Initialize(MapInputName.Data(),debug);
0358   // initialize run number file
0359   if (RunListInputName.CompareTo("")== 0) {
0360       std::cerr << "ERROR: No run list file has been provided, please make sure you do so! Aborting!" << std::endl;
0361       return false;
0362   }
0363   std::map<int,RunInfo> ri=readRunInfosFromFile(RunListInputName.Data(),debug,0);
0364   
0365   // Clean up file headers
0366   TObjArray* tok=ASCIIinputName.Tokenize("/");
0367   // get file name
0368   TString file=((TObjString*)tok->At(tok->GetEntries()-1))->String();
0369   delete tok;
0370   tok=file.Tokenize("_");
0371   TString header=((TObjString*)tok->At(0))->String();
0372   delete tok;
0373   
0374   // Get run number from file & find necessary run infos
0375   TString RunString=header(3,header.Sizeof());
0376   std::map<int,RunInfo>::iterator it=ri.find(RunString.Atoi());
0377   //std::cout<<RunString.Atoi()<<"\t"<<it->second.species<<std::endl;
0378   event.SetRunNumber(RunString.Atoi());
0379   event.SetROtype(ReadOut::Type::Caen);
0380   event.SetBeamName(it->second.species);
0381   event.SetBeamID(it->second.pdg);
0382   event.SetBeamEnergy(it->second.energy);
0383   event.SetVop(it->second.vop);
0384   event.SetVov(it->second.vop-it->second.vbr);
0385   event.SetBeamPosX(it->second.posX);
0386   event.SetBeamPosY(it->second.posY);
0387   calib.SetRunNumber(RunString.Atoi());
0388   calib.SetVop(it->second.vop);
0389   calib.SetVov(it->second.vop-it->second.vbr);  
0390   calib.SetBCCalib(false);
0391   //============================================
0392   // Start decoding file
0393   //============================================
0394   TString aline         = "";
0395   TString versionCAEN   = "";
0396   TObjArray* tokens;
0397   std::map<int,std::vector<Caen> > tmpEvent;
0398   std::map<int,double> tmpTime;
0399   std::map<int,std::vector<Caen> >::iterator itevent;
0400   long tempEvtCounter   = 0;
0401   long writeEvtCounter  = 0;
0402   while(!ASCIIinput.eof()){                                                     // run till end of file is reached and read line by line
0403     aline.ReadLine(ASCIIinput);
0404     if(!ASCIIinput.good()) break;
0405     if(aline[0]=='/'){
0406       if (aline.Contains("File Format Version")){
0407         tokens      = aline.Tokenize(" ");
0408         versionCAEN = ((TObjString*)tokens->At(4))->String();
0409         std::cout << "File version: " << ((TObjString*)tokens->At(4))->String() << std::endl;
0410         
0411         if (!(versionCAEN.CompareTo("3.3") == 0 || versionCAEN.CompareTo("3.1") == 0)){
0412           std::cerr << "This version cannot be converted with the current code, please add the corresponding file format to the converter" << std::endl;
0413           tokens->Clear();
0414           delete tokens;
0415           return false;
0416         }  
0417         
0418         tokens->Clear();
0419         delete tokens;
0420       }
0421       else if(aline.Contains("Run start time")){
0422         tokens    = aline.Tokenize(" ");
0423         int year=((TObjString*)tokens->At(8))->String().Atoi();
0424         int month;
0425         TString Stringmonth=((TObjString*)tokens->At(5))->String();
0426         if(Stringmonth=="Jan") month=1;
0427         else if(Stringmonth=="Feb") month=2;
0428         else if(Stringmonth=="Mar") month=3;
0429         else if(Stringmonth=="Apr") month=4;
0430         else if(Stringmonth=="May") month=5;
0431         else if(Stringmonth=="Jun") month=6;
0432         else if(Stringmonth=="Jul") month=7;
0433         else if(Stringmonth=="Aug") month=8;
0434         else if(Stringmonth=="Sep") month=9;
0435         else if(Stringmonth=="Oct") month=10;
0436         else if(Stringmonth=="Nov") month=11;
0437         else if(Stringmonth=="Dec") month=12;
0438         int day=((TObjString*)tokens->At(6))->String().Atoi();
0439         int hour=((TString)((TObjString*)tokens->At(7))->String()(0,2)).Atoi();
0440         int min=((TString)((TObjString*)tokens->At(7))->String()(3,5)).Atoi();
0441         int sec=((TString)((TObjString*)tokens->At(7))->String()(6,8)).Atoi();
0442         TTimeStamp t=TTimeStamp(year,month,day,hour,min,sec);
0443         event.SetBeginRunTime(t);
0444         calib.SetBeginRunTime(t);
0445         tokens->Clear();
0446         delete tokens;
0447       }
0448       continue;
0449     }
0450     tokens=aline.Tokenize(" \t");
0451     tokens->SetOwner(true);
0452     
0453     if (versionCAEN.CompareTo("3.3") == 0){
0454       int Nfields=tokens->GetEntries();
0455       // std::cout << Nfields << std::endl;
0456       if(((TObjString*)tokens->At(0))->String()=="Brd") {
0457         tokens->Clear();
0458         delete tokens;
0459         continue;
0460       }
0461       //================================================================================
0462       // data format example
0463       // Brd  Ch       LG       HG        Tstamp_us        TrgID        NHits
0464       // 7  00       51       68     10203578.136            0      64
0465       // 7  01       55       75 
0466       //================================================================================
0467       if(Nfields!=7){
0468         std::cout<<"Unexpected number of fields"<<std::endl;
0469         std::cout<<"Expected 7, got: "<<Nfields<<", exit"<<std::endl;
0470         for(int j=0; j<Nfields;j++) std::cout<<j<<"  "<<((TObjString*)tokens->At(j))->String()<<std::endl;
0471         tokens->Clear();
0472         delete tokens;
0473         return -1;
0474       }
0475       int TriggerID=((TObjString*)tokens->At(5))->String().Atoi();
0476       int NHits=((TObjString*)tokens->At(6))->String().Atoi();
0477       double Time = ((TObjString*)tokens->At(4))->String().Atof();
0478       Caen aTile;
0479       int aBoard=((TObjString*)tokens->At(0))->String().Atoi();
0480       int achannel=((TObjString*)tokens->At(1))->String().Atoi();
0481       aTile.SetE(((TObjString*)tokens->At(3))->String().Atoi());//To Test
0482       aTile.SetADCHigh(((TObjString*)tokens->At(3))->String().Atoi());
0483       aTile.SetADCLow (((TObjString*)tokens->At(2))->String().Atoi());
0484       tokens->Clear();
0485       delete tokens;
0486       aTile.SetCellID(setup->GetCellID(aBoard,achannel));
0487       itevent=tmpEvent.find(TriggerID);
0488 
0489       if(itevent!=tmpEvent.end()){
0490         tmpTime[TriggerID]+=Time;
0491         if (aTile.GetCellID() != -1){
0492           itevent->second.push_back(aTile);
0493         } else {
0494           if(debug ==10) std::cout << "cell " << aBoard << "\t" << achannel << " wasn't active according to mapping file!" << std::endl;
0495         }
0496         for(int ich=1; ich<NHits; ich++){
0497           aline.ReadLine(ASCIIinput);
0498           tokens=aline.Tokenize(" ");
0499           tokens->SetOwner(true);
0500           Nfields=tokens->GetEntries();
0501           if(Nfields!=4){
0502             std::cout<<"Expecting 4 fields but read "<<Nfields<<std::endl;
0503             return -1;
0504           }
0505           achannel=((TObjString*)tokens->At(1))->String().Atoi();
0506           aTile.SetE(((TObjString*)tokens->At(3))->String().Atoi());//To Test
0507           aTile.SetADCHigh(((TObjString*)tokens->At(3))->String().Atoi());
0508           aTile.SetADCLow (((TObjString*)tokens->At(2))->String().Atoi());
0509           aTile.SetCellID(setup->GetCellID(aBoard,achannel));
0510 
0511           if (aTile.GetCellID() != -1){
0512             itevent->second.push_back(aTile);
0513           } else {
0514             if(debug ==10) std::cout << "cell " << aBoard << "\t" << achannel << " wasn't active according to mapping file!" << std::endl;
0515           }
0516           tokens->Clear();
0517           delete tokens;
0518         }
0519         // std::cout << itevent->second.size() << "\t" << setup->GetTotalNbChannels() << std::endl;
0520 
0521         if((int)itevent->second.size()==setup->GetTotalNbChannels()/*8*64*/){
0522           //Fill the tree the event is complete and erase from the map
0523           event.SetTimeStamp(tmpTime[TriggerID]/setup->GetNMaxROUnit());
0524           event.SetEventID(itevent->first);
0525           for(std::vector<Caen>::iterator itv=itevent->second.begin(); itv!=itevent->second.end(); ++itv){
0526             event.AddTile(new Caen(*itv));
0527           }
0528           TdataOut->Fill();
0529           writeEvtCounter++;
0530           tmpEvent.erase(itevent);
0531           tmpTime.erase(TriggerID);
0532         } 
0533       } else {
0534         //This is a new event;
0535         tempEvtCounter++;                                                                   // in crease event counts for monitoring of progress
0536         if (tempEvtCounter%5000 == 0 && debug > 0) std::cout << "Converted " <<  tempEvtCounter << " events" << std::endl;
0537         // if (tempEvtCounter > 1000) continue;
0538         std::vector<Caen> vCaen;
0539         if (aTile.GetCellID() != -1){
0540           vCaen.push_back(aTile);
0541         } else {
0542           if(debug ==10) std::cout << "cell " << aBoard << "\t" << achannel << " wasn't active according to mapping file!" << std::endl;
0543         }
0544         for(int ich=1; ich<NHits; ich++){
0545           aline.ReadLine(ASCIIinput);
0546           tokens=aline.Tokenize(" ");
0547           tokens->SetOwner(true);
0548           Nfields=tokens->GetEntries();
0549           if(Nfields!=4){
0550             std::cout<<"Expecting 4 fields but read "<<Nfields<<std::endl;
0551             return -1;
0552           }
0553           achannel=((TObjString*)tokens->At(1))->String().Atoi();
0554           aTile.SetE(((TObjString*)tokens->At(3))->String().Atoi());//To Test
0555           aTile.SetADCHigh(((TObjString*)tokens->At(3))->String().Atoi());
0556           aTile.SetADCLow (((TObjString*)tokens->At(2))->String().Atoi());
0557           aTile.SetCellID(setup->GetCellID(aBoard,achannel));
0558           if (aTile.GetCellID() != -1){
0559             vCaen.push_back(aTile);
0560           } else {
0561             if(debug ==10) std::cout << "cell " << aBoard << "\t" << achannel << " wasn't active according to mapping file!" << std::endl;
0562           }
0563           tokens->Clear();
0564           delete tokens;
0565         }       
0566         tmpTime[TriggerID]=Time;
0567         tmpEvent[TriggerID]=vCaen;
0568         
0569         // std::cout << vCaen.size() << "\t" << setup->GetTotalNbChannels() << std::endl;
0570         if((int)vCaen.size()==setup->GetTotalNbChannels()/*8*64*/){
0571           itevent=tmpEvent.find(TriggerID);
0572           //Fill the tree the event is complete and erase from the map
0573           event.SetTimeStamp(tmpTime[TriggerID]/setup->GetNMaxROUnit());
0574           event.SetEventID(itevent->first);
0575           for(std::vector<Caen>::iterator itv=vCaen.begin(); itv!=vCaen.end(); ++itv){
0576             event.AddTile(new Caen(*itv));
0577           }
0578           TdataOut->Fill();
0579           writeEvtCounter++;
0580           tmpEvent.erase(itevent);
0581           tmpTime.erase(TriggerID);
0582         } 
0583       }
0584     } else if (versionCAEN.CompareTo("3.1") == 0){
0585       int Nfields=tokens->GetEntries();
0586       if(((TObjString*)tokens->At(0))->String()=="Tstamp_us") {
0587         tokens->Clear();
0588         delete tokens;
0589         continue;
0590       }
0591       
0592       //================================================================================
0593       // data format example
0594       //   Tstamp_us        TrgID   Brd  Ch       LG       HG
0595       // 4970514.360            0    01  00       49       46 
0596       //                             01  01       49       35 
0597       //================================================================================
0598       
0599       if(Nfields!=6){
0600         std::cout<<"Unexpected number of fields"<<std::endl;
0601         std::cout<<"Expected 6, got: "<<Nfields<<", exit"<<std::endl;
0602         for(int j=0; j<Nfields;j++) std::cout<<j<<"  "<<((TObjString*)tokens->At(j))->String()<<std::endl;
0603         tokens->Clear();
0604         delete tokens;
0605         return -1;
0606       }
0607       
0608       // std::cout << aline.Data() << std::endl;
0609       int TriggerID = ((TObjString*)tokens->At(1))->String().Atoi();
0610       int NHits     = 64;                                           // temporary fix
0611       double Time   = ((TObjString*)tokens->At(0))->String().Atof();
0612       Caen aTile;
0613       int aBoard    =((TObjString*)tokens->At(2))->String().Atoi();
0614       int achannel  =((TObjString*)tokens->At(3))->String().Atoi();
0615       aTile.SetE(((TObjString*)tokens->At(5))->String().Atoi());//To Test
0616       aTile.SetADCHigh(((TObjString*)tokens->At(5))->String().Atoi());
0617       aTile.SetADCLow (((TObjString*)tokens->At(4))->String().Atoi());
0618       tokens->Clear();
0619       delete tokens;
0620       aTile.SetCellID(setup->GetCellID(aBoard,achannel));
0621       itevent=tmpEvent.find(TriggerID);
0622       
0623       if(itevent!=tmpEvent.end()){
0624         tmpTime[TriggerID]+=Time;
0625         if (aTile.GetCellID() != -1){
0626           itevent->second.push_back(aTile);
0627         } else {
0628           if(debug ==10) std::cout << "cell " << aBoard << "\t" << achannel << " wasn't active according to mapping file!" << std::endl;
0629         }
0630         for(int ich=1; ich<NHits; ich++){
0631             aline.ReadLine(ASCIIinput);
0632             // std::cout << aline.Data() << std::endl;
0633             tokens=aline.Tokenize(" ");
0634             tokens->SetOwner(true);
0635             Nfields=tokens->GetEntries();
0636             
0637             if(Nfields!=4){
0638               std::cout<< "Current line :" << aline.Data() << std::endl;
0639               std::cout<<"Expecting 4 fields but read "<<Nfields<<std::endl;
0640               return -1;
0641             }
0642             achannel=((TObjString*)tokens->At(1))->String().Atoi();
0643             aTile.SetE(((TObjString*)tokens->At(3))->String().Atoi());//To Test
0644             aTile.SetADCHigh(((TObjString*)tokens->At(3))->String().Atoi());
0645             aTile.SetADCLow (((TObjString*)tokens->At(2))->String().Atoi());
0646             aTile.SetCellID(setup->GetCellID(aBoard,achannel));
0647             
0648             if (aTile.GetCellID() != -1){
0649               itevent->second.push_back(aTile);
0650             } else {
0651               if(debug ==10) std::cout << "cell " << aBoard << "\t" << achannel << " wasn't active according to mapping file!" << std::endl;
0652             }
0653             tokens->Clear();
0654             delete tokens;
0655         }
0656         if((int)itevent->second.size()==setup->GetTotalNbChannels()){
0657           //Fill the tree the event is complete and erase from the map
0658           event.SetTimeStamp(tmpTime[TriggerID]/setup->GetNMaxROUnit());
0659           event.SetEventID(itevent->first);
0660           for(std::vector<Caen>::iterator itv=itevent->second.begin(); itv!=itevent->second.end(); ++itv){
0661             event.AddTile(new Caen(*itv));
0662           }
0663           TdataOut->Fill();
0664           writeEvtCounter++;
0665           tmpEvent.erase(itevent);
0666           tmpTime.erase(TriggerID);
0667         } else {
0668           std::cout << "didn't fill" << (int)itevent->second.size()  <<  setup->GetTotalNbChannels()<< std::endl; 
0669         }
0670       } else {
0671         //This is a new event;
0672         tempEvtCounter++;                                                                   // in crease event counts for monitoring of progress
0673         if (tempEvtCounter%5000 == 0 && debug > 0) std::cout << "Converted " <<  tempEvtCounter << " events" << std::endl;
0674         std::vector<Caen> vCaen;
0675         
0676         if (aTile.GetCellID() != -1){
0677           vCaen.push_back(aTile);
0678         } else {
0679           if(debug ==10) std::cout << "cell " << aBoard << "\t" << achannel << " wasn't active according to mapping file!" << std::endl;
0680         }
0681         for(int ich=1; ich<NHits; ich++){
0682           aline.ReadLine(ASCIIinput);
0683           // std::cout << aline.Data() << std::endl;
0684           tokens=aline.Tokenize(" ");
0685           tokens->SetOwner(true);
0686           Nfields=tokens->GetEntries();
0687           if(Nfields!=4){
0688             std::cout<< "Current line :" << aline.Data() << std::endl;
0689             std::cout<<"Expecting 4 fields but read "<<Nfields<<std::endl;
0690             return -1;
0691           }
0692           achannel=((TObjString*)tokens->At(1))->String().Atoi();
0693           aTile.SetE(((TObjString*)tokens->At(3))->String().Atoi());//To Test
0694           aTile.SetADCHigh(((TObjString*)tokens->At(3))->String().Atoi());
0695           aTile.SetADCLow (((TObjString*)tokens->At(2))->String().Atoi());
0696           aTile.SetCellID(setup->GetCellID(aBoard,achannel));
0697           if (aTile.GetCellID() != -1){
0698             vCaen.push_back(aTile);
0699           } else {
0700             if(debug ==10) std::cout << "cell " << aBoard << "\t" << achannel << " wasn't active according to mapping file!" << std::endl;
0701           }
0702           tokens->Clear();
0703           delete tokens;
0704         }
0705         tmpTime[TriggerID]=Time;
0706         tmpEvent[TriggerID]=vCaen;
0707         
0708         if((int)vCaen.size()==setup->GetTotalNbChannels()/*8*64*/){
0709           itevent=tmpEvent.find(TriggerID);
0710           //Fill the tree the event is complete and erase from the map
0711           event.SetTimeStamp(tmpTime[TriggerID]/setup->GetNMaxROUnit());
0712           event.SetEventID(itevent->first);
0713           for(std::vector<Caen>::iterator itv=vCaen.begin(); itv!=vCaen.end(); ++itv){
0714             event.AddTile(new Caen(*itv));
0715           }
0716           TdataOut->Fill();
0717           writeEvtCounter++;
0718           tmpEvent.erase(itevent);
0719           tmpTime.erase(TriggerID);
0720         }
0721         
0722       }      
0723     }
0724   } // finished reading in file
0725   // 
0726   if (debug > 0) std::cout << "Converted a total of " <<  tempEvtCounter << " events" << std::endl;
0727   
0728   //============================================
0729   // Fill & write all trees to output file 
0730   //============================================
0731   RootOutput->cd();
0732   // setup 
0733   RootSetupWrapper rswtmp=RootSetupWrapper(setup);
0734   rsw=rswtmp;
0735   TsetupOut->Fill();
0736   // calib
0737   TcalibOut->Fill();
0738   TcalibOut->Write();
0739   // event data
0740   TdataOut->Fill();
0741   TsetupOut->Write();
0742   TdataOut->Write();
0743 
0744   std::cout << "Events written to file: " << writeEvtCounter << std::endl;
0745   if (writeEvtCounter < 2){
0746     std::cout << "ERROR: Only " << writeEvtCounter << " events were written, something didn't go right, please check your mapping file!" << std::endl; 
0747   }
0748   RootOutput->Close();
0749   return true;
0750 }
0751 
0752 
0753 
0754 // ****************************************************************************
0755 // convert already processed root file from CAEN output to new root format
0756 // ****************************************************************************
0757 bool Analyses::ConvertOldRootFile2Root(void){
0758   //============================================
0759   //Init first
0760   //============================================
0761   // initialize setup
0762   if (MapInputName.CompareTo("")== 0) {
0763       std::cerr << "ERROR: No mapping file has been provided, please make sure you do so! Aborting!" << std::endl;
0764       return false;
0765   }
0766   setup->Initialize(MapInputName.Data(),debug);
0767   // initialize run number file
0768   if (RunListInputName.CompareTo("")== 0) {
0769       std::cerr << "ERROR: No run list file has been provided, please make sure you do so! Aborting!" << std::endl;
0770       return false;
0771   }
0772   std::map<int,RunInfo> ri=readRunInfosFromFile(RunListInputName.Data(),debug,0);
0773   
0774   // Clean up file headers
0775   TObjArray* tok=ASCIIinputName.Tokenize("/");
0776   // get file name
0777   TString file=((TObjString*)tok->At(tok->GetEntries()-1))->String();
0778   delete tok;
0779   tok=file.Tokenize("_");
0780   TString header=((TObjString*)tok->At(0))->String();
0781   delete tok;
0782   
0783   // Get run number from file & find necessary run infos
0784   TString RunString=header(3,header.Sizeof());
0785   std::map<int,RunInfo>::iterator it=ri.find(RunString.Atoi());
0786   //std::cout<<RunString.Atoi()<<"\t"<<it->second.species<<std::endl;
0787   event.SetRunNumber(RunString.Atoi());
0788   event.SetROtype(ReadOut::Type::Caen);
0789   event.SetBeamName(it->second.species);
0790   event.SetBeamID(it->second.pdg);
0791   event.SetBeamEnergy(it->second.energy);
0792   event.SetVop(it->second.vop);
0793   event.SetVov(it->second.vop-it->second.vbr);
0794   event.SetBeamPosX(it->second.posX);
0795   event.SetBeamPosY(it->second.posY);
0796   calib.SetRunNumber(RunString.Atoi());
0797   calib.SetVop(it->second.vop);
0798   calib.SetVov(it->second.vop-it->second.vbr);  
0799   calib.SetBCCalib(false);
0800     // load tree
0801   TChain *const tt_event = new TChain("tree");
0802   if (ASCIIinputName.EndsWith(".root")) {                     // are we loading a single root tree?
0803       std::cout << "loading a single root file" << std::endl;
0804       tt_event->AddFile(ASCIIinputName);
0805       TFile testFile(ASCIIinputName.Data());
0806       if (testFile.IsZombie()){
0807         std::cout << Form("The file %s is not a root file or doesn't exit, please fix the file path", ASCIIinputName.Data()) << std::endl;
0808         return false;
0809       }
0810 
0811   } else {
0812       std::cout << "please try again this isn't a root file" << std::endl;
0813   } 
0814   if(!tt_event){ std::cout << "tree not found... returning!"<< std::endl; return false;}
0815 
0816   // Define tree variables
0817   Long64_t gTrID;
0818   Double_t gTRtimeStamp;
0819   const int gMaxChannels = 64;
0820   Long64_t* gBoard          = new Long64_t[gMaxChannels];
0821   Long64_t* gChannel        = new Long64_t[gMaxChannels];
0822   Long64_t* gLG             = new Long64_t[gMaxChannels];
0823   Long64_t* gHG             = new Long64_t[gMaxChannels];
0824 
0825   if (tt_event->GetBranchStatus("t_stamp") ){
0826     tt_event->SetBranchAddress("trgid",            &gTrID);
0827     tt_event->SetBranchAddress("t_stamp",          &gTRtimeStamp);
0828     tt_event->SetBranchAddress("board",            gBoard);
0829     tt_event->SetBranchAddress("channel",          gChannel);
0830     tt_event->SetBranchAddress("LG",               gLG);
0831     tt_event->SetBranchAddress("HG",               gHG);
0832   }
0833   
0834   Long64_t nEntriesTree                 = tt_event->GetEntries();
0835   std::cout << "Number of events in tree: " << nEntriesTree << std::endl;
0836 
0837   std::map<int,std::vector<Caen> > tmpEvent;
0838   std::map<int,double> tmpTime;
0839   for (Long64_t i=0; i<nEntriesTree;i++) {
0840     // load current event
0841     tt_event->GetEntry(i);  
0842     if (i%5000 == 0 && debug > 0) std::cout << "Converted " <<  i << " events" << std::endl;    
0843     int TriggerID = gTrID;
0844     double Time   = gTRtimeStamp;
0845     std::vector<Caen> vCaen;
0846     
0847     for(Int_t ch=0; ch<gMaxChannels; ch++){   
0848       Caen aTile;
0849       int aBoard      = gBoard[ch];
0850       int achannel    = gChannel[ch];
0851       aTile.SetE(gHG[ch]);//To Test
0852       aTile.SetADCHigh(gHG[ch]);
0853       aTile.SetADCLow (gLG[ch]);
0854       aTile.SetCellID(setup->GetCellID(aBoard,achannel));
0855       if (aTile.GetCellID() != -1){
0856         vCaen.push_back(aTile);
0857       } else {
0858         if(debug ==10) std::cout << "cell " << aBoard << "\t" << achannel << " wasn't active according to mapping file!" << std::endl;
0859       }
0860     }
0861     
0862      if((int)vCaen.size()==setup->GetTotalNbChannels()){
0863       //Fill the tree the event is complete and erase from the map
0864       event.SetTimeStamp(Time);
0865       event.SetEventID(TriggerID);
0866       for(std::vector<Caen>::iterator itv=vCaen.begin(); itv!=vCaen.end(); ++itv){
0867         event.AddTile(new Caen(*itv));
0868       }
0869       TdataOut->Fill();
0870       vCaen.clear();
0871     }
0872   } // finished reading the file
0873 
0874   // 
0875   if (debug > 0) std::cout << "Converted a total of " <<  nEntriesTree << " events" << std::endl;
0876   
0877   //============================================
0878   // Fill & write all trees to output file 
0879   //============================================
0880   RootOutput->cd();
0881   // setup 
0882   RootSetupWrapper rswtmp=RootSetupWrapper(setup);
0883   rsw=rswtmp;
0884   TsetupOut->Fill();
0885   // calib
0886   TcalibOut->Fill();
0887   TcalibOut->Write();
0888   // event data
0889   TdataOut->Fill();
0890   TsetupOut->Write();
0891   TdataOut->Write();
0892 
0893   RootOutput->Close();
0894   return true;
0895 }
0896 
0897 
0898 
0899 // ****************************************************************************
0900 // extract pedestral from dedicated data run, no other data in run 
0901 // ****************************************************************************
0902 bool Analyses::GetPedestal(void){
0903   std::cout<<"GetPedestal for maximium "<< setup->GetMaxCellID() << " cells" <<std::endl;
0904   
0905   std::map<int,RunInfo> ri=readRunInfosFromFile(RunListInputName.Data(),debug,0);
0906   
0907   int maxChannelPerLayer = (setup->GetNMaxColumn()+1)*(setup->GetNMaxRow()+1);
0908   ROOT::Math::MinimizerOptions::SetDefaultMinimizer("Minuit2", "Migrad");  
0909   
0910   // create HG and LG histo's per channel
0911   TH2D* hspectraHGvsCellID      = new TH2D( "hspectraHG_vsCellID","ADC spectrum High Gain vs CellID; cell ID; ADC_{HG} (arb. units) ",
0912                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5, 4000,0,4000);
0913   hspectraHGvsCellID->SetDirectory(0);
0914   TH2D* hspectraLGvsCellID      = new TH2D( "hspectraLG_vsCellID","ADC spectrum Low Gain vs CellID; cell ID; ADC_{LG} (arb. units)  ",
0915                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5, 4000,0,4000);
0916   hspectraLGvsCellID->SetDirectory(0);
0917   TH1D* hMeanPedHGvsCellID      = new TH1D( "hMeanPedHG_vsCellID","mean Ped High Gain vs CellID ; cell ID; #mu_{noise, HG} (arb. units) ",
0918                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5);
0919   hMeanPedHGvsCellID->SetDirectory(0);
0920   TH2D* hspectraHGMeanVsLayer   = new TH2D( "hspectraHGMeanVsLayer","Mean Ped High Gain vs CellID; layer; brd channel; #mu_{Ped HG} (arb. units) ",
0921                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
0922   hspectraHGMeanVsLayer->SetDirectory(0);
0923   TH2D* hspectraHGSigmaVsLayer  = new TH2D( "hspectraHGSigmaVsLayer","Mean Ped High Gain vs CellID; layer; brd channel; #sigma_{Ped HG} (arb. units) ",
0924                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
0925   hspectraHGSigmaVsLayer->SetDirectory(0);
0926   TH1D* hMeanPedLGvsCellID      = new TH1D( "hMeanPedLG_vsCellID","mean Ped Low Gain vs CellID ; cell ID; #mu_{noise, LG} (arb. units) ",
0927                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5);
0928   hMeanPedLGvsCellID->SetDirectory(0);
0929   TH2D* hspectraLGMeanVsLayer   = new TH2D( "hspectraLGMeanVsLayer","Mean Ped Low Gain vs CellID; layer; brd channel; #mu_{PED LG} (arb. units) ",
0930                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
0931   hspectraLGMeanVsLayer->SetDirectory(0);
0932   TH2D* hspectraLGSigmaVsLayer  = new TH2D( "hspectraLGSigmaVsLayer","Mean Ped Low Gain vs CellID; layer; brd channel; #sigma_{Ped LG} (arb. units)",
0933                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
0934   hspectraLGSigmaVsLayer->SetDirectory(0);
0935 
0936   TH2D* hPedMeanHGvsLG          = new TH2D( "hPedMeanHGvsLG","Mean Ped High Gain vs Low Gain; #mu_{noise, HG} (arb. units); #mu_{noise, LG} (arb. units)", 500, 0, 250, 500, 0, 250);
0937   hPedMeanHGvsLG->SetDirectory(0);
0938   
0939   std::map<int,TileSpectra> hSpectra;
0940   std::map<int, TileSpectra>::iterator ithSpectra;
0941 
0942   std::cout << "Additional Output with histos being created: " << RootOutputNameHist.Data() << std::endl;
0943   if(Overwrite){
0944     std::cout << "recreating file with hists" << std::endl;
0945     RootOutputHist = new TFile(RootOutputNameHist.Data(),"RECREATE");
0946   } else{
0947     std::cout << "newly creating file with hists" << std::endl;
0948     RootOutputHist = new TFile(RootOutputNameHist.Data(),"CREATE");
0949   }
0950   // entering histoOutput file
0951   RootOutputHist->mkdir("IndividualCells");
0952   RootOutputHist->cd("IndividualCells");
0953 
0954   RootOutput->cd();
0955   // Event loop to fill histograms & output tree
0956   std::cout << "N max layers: " << setup->GetNMaxLayer() << "\t columns: " <<  setup->GetNMaxColumn() << "\t row: " << setup->GetNMaxRow() << "\t module: " <<  setup->GetNMaxModule() << std::endl;
0957   if(TcalibIn) TcalibIn->GetEntry(0);
0958   // check whether calib should be overwritten based on external text file
0959   if (OverWriteCalib){
0960     calib.ReadCalibFromTextFile(ExternalCalibFile,debug);
0961   }
0962   
0963   int evts=TdataIn->GetEntries();
0964   int runNr = -1;
0965   if (maxEvents == -1){
0966     maxEvents = evts;
0967   } else {
0968     std::cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl;
0969     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;
0970     std::cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl;
0971   }
0972   ReadOut::Type typeRO = ReadOut::Type::Caen;
0973   
0974   for(int i=0; i<evts && i < maxEvents; i++){
0975     TdataIn->GetEntry(i);
0976     if (i == 0){
0977       runNr   = event.GetRunNumber();
0978       typeRO  = event.GetROtype();
0979       std::cout<< "original run numbers calib: "<<calib.GetRunNumber() << "\t" << calib.GetRunNumberPed() << "\t" << calib.GetRunNumberMip() << std::endl;
0980       if (typeRO == ReadOut::Type::Caen) std::cout << "Read-out type CAEN" << std::endl;
0981       else  std::cout << "Read-out type HGCROC" << std::endl;
0982       calib.SetRunNumberPed(runNr);
0983       calib.SetBeginRunTimePed(event.GetBeginRunTimeAlt());
0984       std::cout<< "reset run numbers calib: "<< calib.GetRunNumber() << "\t" << calib.GetRunNumberPed() << "\t" << calib.GetRunNumberMip() << std::endl;
0985     }
0986     if (i%5000 == 0&& i > 0 && debug > 0) std::cout << "Reading " <<  i << "/" << evts << " events"<< std::endl;
0987     for(int j=0; j<event.GetNTiles(); j++){
0988       if (typeRO == ReadOut::Type::Caen) {
0989         Caen* aTile=(Caen*)event.GetTile(j);
0990         if (i == 0 && debug > 2){
0991           std::cout << ((TString)setup->DecodeCellID(aTile->GetCellID())).Data() << std::endl;
0992         }
0993         ithSpectra=hSpectra.find(aTile->GetCellID());
0994         if(ithSpectra!=hSpectra.end()){
0995           ithSpectra->second.Fill(aTile->GetADCLow(),aTile->GetADCHigh());
0996         } else {
0997           RootOutputHist->cd("IndividualCells");
0998           hSpectra[aTile->GetCellID()]=TileSpectra("1stExtraction",aTile->GetCellID(),calib.GetTileCalib(aTile->GetCellID()), event.GetROtype(), debug);
0999           hSpectra[aTile->GetCellID()].Fill(aTile->GetADCLow(),aTile->GetADCHigh());
1000           RootOutput->cd();
1001         }
1002         hspectraHGvsCellID->Fill(aTile->GetCellID(), aTile->GetADCHigh());
1003         hspectraLGvsCellID->Fill(aTile->GetCellID(), aTile->GetADCLow());
1004       }
1005       else if (typeRO == ReadOut::Type::Hgcroc){ // Process HGCROC Data
1006         Hgcroc* aTile=(Hgcroc*)event.GetTile(j);
1007         if (i == 0 && debug > 2){
1008           std::cout << ((TString)setup->DecodeCellID(aTile->GetCellID())).Data() << std::endl;
1009         }
1010         ithSpectra=hSpectra.find(aTile->GetCellID());
1011         if(ithSpectra!=hSpectra.end()){
1012           // Get the tile spectra if it already exists
1013           ithSpectra->second.FillExt(aTile->GetPedestal(),aTile->GetPedestal(), 0., 0.);
1014           ithSpectra->second.FillWaveform(aTile->GetADCWaveform());          
1015         } else {
1016           // Make a new tile spectra if it isn't found
1017           RootOutputHist->cd("IndividualCells");
1018           hSpectra[aTile->GetCellID()]= TileSpectra("1stExtraction", 2, aTile->GetCellID(), calib.GetTileCalib(aTile->GetCellID()), event.GetROtype(), debug);
1019           
1020           hSpectra[aTile->GetCellID()].FillExt(aTile->GetPedestal(),aTile->GetPedestal(), 0., 0.);
1021           hSpectra[aTile->GetCellID()].FillWaveform(aTile->GetADCWaveform());
1022 
1023           RootOutput->cd();
1024         }
1025         // std::cout << "Cell ID: " << aTile->GetCellID() << ", Tile E: " << aTile->GetPedestal() << std::endl;
1026         hspectraHGvsCellID->Fill(aTile->GetCellID(), aTile->GetPedestal());
1027         hspectraLGvsCellID->Fill(aTile->GetCellID(), aTile->GetPedestal());
1028       }
1029     }
1030     RootOutput->cd();
1031     TdataOut->Fill();
1032   }
1033   
1034   // fit pedestal
1035   double* parameters=new double[8];
1036   bool isGood  = true;
1037   bool isGood2 = true;
1038 
1039   for(ithSpectra=hSpectra.begin(); ithSpectra!=hSpectra.end(); ++ithSpectra){
1040     if ( debug > 2) std::cout << ((TString)setup->DecodeCellID(ithSpectra->second.GetCellID())).Data() << std::endl;
1041     // std::cerr << "Fitting noise for cell ID: " << ithSpectra->second.GetCellID() << std::endl;
1042     isGood=ithSpectra->second.FitNoise(parameters, yearData, false);
1043     if (!isGood && !(typeRO == ReadOut::Type::Hgcroc)) {
1044       std::cerr << "Noise fit failed for cell ID: " << ithSpectra->second.GetCellID() << std::endl;
1045       continue;
1046     }
1047 
1048     if (typeRO == ReadOut::Type::Hgcroc){
1049       isGood2=ithSpectra->second.FitPedConstWage(debug);
1050 
1051       if (!isGood && !isGood2) {
1052         std::cerr << "both noise fits failed " << ithSpectra->second.GetCellID() << std::endl;
1053         continue;
1054       } else {
1055         if (!isGood2) std::cerr << "Noise-wave form fit failed for cell ID: " << ithSpectra->second.GetCellID() << std::endl;
1056         if (!isGood){
1057           std::cerr << "1D Noise fit failed for cell ID: " << ithSpectra->second.GetCellID() << std::endl;
1058           parameters[4] = -1;
1059           parameters[5] = -1;
1060           parameters[6] = -1;
1061           parameters[7] = -1;
1062         }
1063       }
1064     }
1065     
1066     
1067     int layer     = setup->GetLayer(ithSpectra->second.GetCellID());
1068     int chInLayer = setup->GetChannelInLayer(ithSpectra->second.GetCellID());
1069 
1070     hMeanPedHGvsCellID->SetBinContent(hMeanPedHGvsCellID->FindBin(ithSpectra->second.GetCellID()), parameters[4]);
1071     hMeanPedHGvsCellID->SetBinError  (hMeanPedHGvsCellID->FindBin(ithSpectra->second.GetCellID()), parameters[6]);
1072     hspectraHGMeanVsLayer->SetBinContent(hspectraHGMeanVsLayer->FindBin(layer,chInLayer), parameters[4]);
1073     hspectraHGMeanVsLayer->SetBinError(hspectraHGMeanVsLayer->FindBin(layer,chInLayer), parameters[5]);
1074     hspectraHGSigmaVsLayer->SetBinContent(hspectraHGSigmaVsLayer->FindBin(layer,chInLayer), parameters[6]);
1075     hspectraHGSigmaVsLayer->SetBinError(hspectraHGSigmaVsLayer->FindBin(layer,chInLayer), parameters[7]);
1076     if (typeRO == ReadOut::Type::Caen){
1077       hMeanPedLGvsCellID->SetBinContent(hMeanPedLGvsCellID->FindBin(ithSpectra->second.GetCellID()), parameters[0]);
1078       hMeanPedLGvsCellID->SetBinError  (hMeanPedLGvsCellID->FindBin(ithSpectra->second.GetCellID()), parameters[2]);
1079       hspectraLGMeanVsLayer->SetBinContent(hspectraLGMeanVsLayer->FindBin(layer,chInLayer), parameters[0]);
1080       hspectraLGMeanVsLayer->SetBinError(hspectraLGMeanVsLayer->FindBin(layer,chInLayer), parameters[1]);
1081       hspectraLGSigmaVsLayer->SetBinContent(hspectraLGSigmaVsLayer->FindBin(layer,chInLayer), parameters[2]);
1082       hspectraLGSigmaVsLayer->SetBinError(hspectraLGSigmaVsLayer->FindBin(layer,chInLayer), parameters[3]);
1083       
1084       hPedMeanHGvsLG->Fill(parameters[4],parameters[0]);
1085     } else if (isGood2 && typeRO == ReadOut::Type::Hgcroc){
1086       hMeanPedLGvsCellID->SetBinContent(hMeanPedLGvsCellID->FindBin(ithSpectra->second.GetCellID()), calib.GetPedestalMeanL(ithSpectra->second.GetCellID()));
1087       hMeanPedLGvsCellID->SetBinError  (hMeanPedLGvsCellID->FindBin(ithSpectra->second.GetCellID()), calib.GetPedestalSigL(ithSpectra->second.GetCellID()));
1088       hspectraLGMeanVsLayer->SetBinContent(hspectraLGMeanVsLayer->FindBin(layer,chInLayer), calib.GetPedestalMeanL(ithSpectra->second.GetCellID()));
1089       hspectraLGMeanVsLayer->SetBinError(hspectraLGMeanVsLayer->FindBin(layer,chInLayer), calib.GetPedestalSigL(ithSpectra->second.GetCellID()));
1090     }
1091   }
1092   
1093   RootOutput->cd();
1094   // write output tree (copy of raw data)
1095   TdataOut->Write();
1096   // write setup tree (copy of raw data)
1097   TsetupIn->CloneTree()->Write();
1098   
1099   if (IsCalibSaveToFile()){
1100     TString fileCalibPrint = RootOutputName;
1101     fileCalibPrint         = fileCalibPrint.ReplaceAll(".root","_calib.txt");
1102     calib.PrintCalibToFile(fileCalibPrint);
1103   }
1104   
1105   TcalibOut->Fill();
1106   TcalibOut->Write();
1107   RootOutput->Write();
1108   RootOutput->Close();
1109   
1110   // entering histoOutput file
1111   //RootOutputHist->cd();
1112   //RootOutputHist->mkdir("IndividualCells");
1113   RootOutputHist->cd("IndividualCells");
1114   for(ithSpectra=hSpectra.begin(); ithSpectra!=hSpectra.end(); ++ithSpectra){
1115     ithSpectra->second.Write(true);
1116   }
1117   RootOutputHist->cd();
1118   hspectraHGvsCellID->Write();
1119   hMeanPedHGvsCellID->Write();
1120   hspectraHGMeanVsLayer->Write();
1121   hspectraHGSigmaVsLayer->Write();
1122 
1123   if (typeRO == ReadOut::Type::Caen){
1124     hspectraLGvsCellID->Write();
1125     hMeanPedLGvsCellID->Write();
1126     hspectraLGMeanVsLayer->Write();
1127     hspectraLGSigmaVsLayer->Write();
1128     hPedMeanHGvsLG->Write();
1129   } else {
1130     hMeanPedLGvsCellID->GetYaxis()->SetTitle("#mu_{Ped ADC, wave} (arb. units)");
1131     hMeanPedLGvsCellID->Write("hMeanPedWave_vsCellID");
1132     hspectraLGMeanVsLayer->GetZaxis()->SetTitle("#mu_{Ped ADC, wave} (arb. units)");
1133     hspectraLGMeanVsLayer->Write("hspectraPedWaveMeanVsLayer");
1134   }
1135   // fill calib tree & write it
1136   // close open root files
1137   RootOutputHist->Write();
1138   RootOutputHist->Close();
1139 
1140   RootInput->Close();
1141   
1142   
1143   // Get run info object
1144   std::map<int,RunInfo>::iterator it=ri.find(runNr);
1145   
1146   // create directory for plot output
1147   TString outputDirPlots = GetPlotOutputDir();
1148   gSystem->Exec("mkdir -p "+outputDirPlots);
1149   
1150   double averagePedMeanHG = calib.GetAveragePedestalMeanHigh();
1151   double averagePedSigHG  = calib.GetAveragePedestalSigHigh();
1152   double averagePedMeanLG = calib.GetAveragePedestalMeanLow();
1153   double averagePedSigLG  = calib.GetAveragePedestalSigLow();
1154 
1155   //**********************************************************************
1156   // Create canvases for channel overview plotting
1157   //**********************************************************************
1158   Double_t textSizeRel = 0.035;
1159   StyleSettingsBasics("pdf");
1160   SetPlotStyle();
1161   
1162   TCanvas* canvas2DCorr = new TCanvas("canvasCorrPlots","",0,0,1450,1200);  // gives the page size
1163   DefaultCancasSettings( canvas2DCorr, 0.08, 0.13, 0.045, 0.07);
1164   canvas2DCorr->SetLogz();
1165   
1166   if (typeRO == ReadOut::Type::Hgcroc) hspectraHGvsCellID->GetYaxis()->SetTitle("Pedestal ADC (arb units)");
1167   PlotSimple2D( canvas2DCorr, hspectraHGvsCellID, 300, setup->GetMaxCellID()+1, textSizeRel, Form("%s/HG_Noise.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 5, kFALSE, "colz", true);
1168   
1169   
1170   if (typeRO == ReadOut::Type::Caen){
1171     PlotSimple2D( canvas2DCorr, hspectraLGvsCellID, 300, setup->GetMaxCellID()+1, textSizeRel, Form("%s/LG_Noise.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 5, kFALSE, "colz", true);
1172   } 
1173   
1174   canvas2DCorr->SetLogz(0);
1175   PlotSimple2D( canvas2DCorr, hspectraHGMeanVsLayer, -10000, -10000, textSizeRel, Form("%s/HG_NoiseMean.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 5, kFALSE, "colz", true, Form("#LT#mu_{HG}#GT = %2.2f", averagePedMeanHG));
1176   PlotSimple2D( canvas2DCorr, hspectraHGSigmaVsLayer,-10000, -10000, textSizeRel, Form("%s/HG_NoiseSigma.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 5,  kFALSE, "colz", true, Form("#LT#sigma_{HG}#GT = %2.2f", averagePedSigHG));
1177 
1178   if (typeRO == ReadOut::Type::Caen){
1179     PlotSimple2D( canvas2DCorr, hspectraLGMeanVsLayer, -10000, -10000, textSizeRel, Form("%s/LG_NoiseMean.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 5, kFALSE, "colz", true, Form("#LT#mu_{LG}#GT = %2.2f", averagePedMeanLG));
1180     PlotSimple2D( canvas2DCorr, hspectraLGSigmaVsLayer, -10000, -10000, textSizeRel, Form("%s/LG_NoiseSigma.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 5, kFALSE, "colz", true, Form("#LT#sigma_{LG}#GT = %2.2f", averagePedSigLG));
1181   
1182     PlotSimple2D( canvas2DCorr, hPedMeanHGvsLG, -10000, -10000, textSizeRel, Form("%s/PedMean_HG_LG.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 5, kFALSE, "colz", true, "");
1183   } else {
1184     PlotSimple2D( canvas2DCorr, hspectraLGMeanVsLayer, -10000, -10000, textSizeRel, Form("%s/PedWave_NoiseMean.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 5, kFALSE, "colz", true, Form("#LT#mu_{wave}#GT = %2.2f", averagePedMeanLG));
1185   }
1186   //***********************************************************************************************************
1187   //********************************* 8 Panel overview plot  **************************************************
1188   //***********************************************************************************************************
1189   //*****************************************************************
1190     // Test beam geometry (beam coming from viewer)
1191     //===========================================================
1192     //||    8 (4)    ||    7 (5)   ||    6 (6)   ||    5 (7)   ||  row 0
1193     //===========================================================
1194     //||    1 (0)    ||    2 (1)   ||    3 (2)   ||    4 (3)   ||  row 1
1195     //===========================================================
1196     //    col 0     col 1       col 2     col  3
1197     // rebuild pad geom in similar way (numbering -1)
1198   //*****************************************************************
1199   TCanvas* canvas8Panel;
1200   TPad* pad8Panel[8];
1201   Double_t topRCornerX[8];
1202   Double_t topRCornerY[8];
1203   Int_t textSizePixel = 30;
1204   Double_t relSize8P[8];
1205   CreateCanvasAndPadsFor8PannelTBPlot(canvas8Panel, pad8Panel,  topRCornerX, topRCornerY, relSize8P, textSizePixel);
1206 
1207   TCanvas* canvas8PanelProf;
1208   TPad* pad8PanelProf[8];
1209   Double_t topRCornerXProf[8];
1210   Double_t topRCornerYProf[8];
1211   Double_t relSize8PProf[8];
1212   CreateCanvasAndPadsFor8PannelTBPlot(canvas8PanelProf, pad8PanelProf,  topRCornerXProf, topRCornerYProf, relSize8PProf, textSizePixel, 0.045, "Prof", false);
1213  
1214   
1215   for (Int_t l = 0; l < setup->GetNMaxLayer()+1; l++){    
1216     PlotNoiseWithFitsFullLayer (canvas8Panel,pad8Panel, topRCornerX, topRCornerY, relSize8P, textSizePixel, 
1217                                 hSpectra, setup, true, 0, 275, 1.2, l, 0,
1218                                 Form("%s/Noise_HG_Layer%02d.%s" ,outputDirPlots.Data(), l, plotSuffix.Data()), it->second);
1219     if (typeRO == ReadOut::Type::Caen){
1220       PlotNoiseWithFitsFullLayer (canvas8Panel,pad8Panel, topRCornerX, topRCornerY, relSize8P, textSizePixel, 
1221                                 hSpectra, setup, false, 0, 275, 1.2, l, 0,
1222                                 Form("%s/Noise_LG_Layer%02d.%s" ,outputDirPlots.Data(), l, plotSuffix.Data()), it->second);
1223     } else if (typeRO == ReadOut::Type::Hgcroc){
1224       PlotCorr2DFullLayer(canvas8PanelProf,pad8PanelProf, topRCornerXProf, topRCornerYProf, relSize8PProf, textSizePixel, hSpectra, 0, it->second.samples+1, 300, l, 0,
1225                                   Form("%s/Waveform_Layer%02d.%s" ,outputDirPlots.Data(), l, plotSuffix.Data()), it->second);
1226  
1227     }
1228   }
1229 
1230   
1231   return true;
1232 }
1233 
1234 // ****************************************************************************
1235 // extract pedestral from dedicated data run, no other data in run 
1236 // ****************************************************************************
1237 bool Analyses::TransferCalib(void){
1238   std::cout<<"Transfer calib"<<std::endl;
1239   int evts=TdataIn->GetEntries();
1240 
1241   std::map<int,RunInfo> ri=readRunInfosFromFile(RunListInputName.Data(),debug,0);
1242   
1243   std::map<int,TileSpectra> hSpectra;
1244   std::map<int, TileSpectra>::iterator ithSpectra;
1245   TcalibIn->GetEntry(0);
1246   // check whether calib should be overwritten based on external text file
1247   if (OverWriteCalib){
1248     calib.ReadCalibFromTextFile(ExternalCalibFile,debug);
1249   }
1250   
1251   int runNr = -1;
1252 
1253   std::map<int,short> bcmap;
1254   std::map<int,short>::iterator bcmapIte;
1255   if (CalcBadChannel == 1)
1256     bcmap = ReadExternalBadChannelMap();
1257 
1258   // *******************************************************************
1259   // ***** create additional output hist *******************************
1260   // *******************************************************************
1261   std::cout << "Additional Output with histos being created: " << RootOutputNameHist.Data() << std::endl;
1262   if(Overwrite){
1263     std::cout << "recreating file with hists" << std::endl;
1264     RootOutputHist = new TFile(RootOutputNameHist.Data(),"RECREATE");
1265   } else{
1266     std::cout << "newly creating file with hists" << std::endl;
1267     RootOutputHist = new TFile(RootOutputNameHist.Data(),"CREATE");
1268   }
1269   TH2D* hspectraADCvsCellID      = new TH2D( "hspectraADCvsCellID","ADC spectrums CellID; cell ID; ADC (arb. units) ",
1270                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5, 1024,0,1024);
1271   hspectraADCvsCellID->SetDirectory(0);
1272   TH2D* hspectraTOTvsCellID      = new TH2D( "hspectraTOTvsCellID","TOT spectrums CellID; cell ID; ADC (arb. units) ",
1273                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5, 4096,0,4096);
1274   hspectraTOTvsCellID->SetDirectory(0);
1275   
1276   TH2D* hspectraADCPedvsCellID      = new TH2D( "hspectraADCPedvsCellID","Pedestal ADC spectrums CellID; cell ID; ADC (arb. units) ",
1277                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5, 1024,0,1024);
1278   hspectraADCPedvsCellID->SetDirectory(0);
1279   
1280   RootOutputHist->mkdir("IndividualCells");
1281   RootOutputHist->cd("IndividualCells");
1282   
1283   
1284   
1285   
1286   ROOT::Math::MinimizerOptions::SetDefaultMinimizer("Minuit2", "Migrad");  
1287   if (maxEvents == -1){
1288     maxEvents = evts;
1289   } else {
1290     std::cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl;
1291     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;
1292     std::cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl;
1293   }
1294   ReadOut::Type typeRO = ReadOut::Type::Caen;
1295   for(int i=0; i<evts && i < maxEvents ; i++){
1296     if (i%5000 == 0&& i > 0 && debug > 0) std::cout << "Reading " <<  i << "/" << evts << " events"<< std::endl;
1297     TdataIn->GetEntry(i);
1298     if (i == 0){
1299       runNr   = event.GetRunNumber();
1300       typeRO  = event.GetROtype();
1301       std::cout<< "original run numbers calib: "<<calib.GetRunNumber() << "\t" << calib.GetRunNumberPed() << "\t" << calib.GetRunNumberMip() << std::endl;
1302       calib.SetRunNumber(runNr);
1303       calib.SetBeginRunTime(event.GetBeginRunTimeAlt());
1304       std::cout<< "reset run numbers calib: "<< calib.GetRunNumber() << "\t" << calib.GetRunNumberPed() << "\t" << calib.GetRunNumberMip() << std::endl;
1305     }
1306    
1307     // if (CalcBadChannel > 0 || ExtPlot > 0){
1308       for(int j=0; j<event.GetNTiles(); j++){
1309         if (typeRO == ReadOut::Type::Caen) {
1310           Caen* aTile=(Caen*)event.GetTile(j);
1311           ithSpectra=hSpectra.find(aTile->GetCellID());
1312           double hgCorr = aTile->GetADCHigh()-calib.GetPedestalMeanH(aTile->GetCellID());
1313           double lgCorr = aTile->GetADCLow()-calib.GetPedestalMeanL(aTile->GetCellID());
1314           
1315           if(ithSpectra!=hSpectra.end()){
1316             ithSpectra->second.FillExt(lgCorr,hgCorr, 0., 0.);
1317           } else {
1318             RootOutputHist->cd("IndividualCells");
1319             hSpectra[aTile->GetCellID()]=TileSpectra("ped",2,aTile->GetCellID(),calib.GetTileCalib(aTile->GetCellID()),event.GetROtype(),debug);
1320             hSpectra[aTile->GetCellID()].FillExt(lgCorr,hgCorr, 0., 0.);
1321             RootOutput->cd();
1322           }
1323         } else if (typeRO == ReadOut::Type::Hgcroc) {
1324           Hgcroc* aTile=(Hgcroc*)event.GetTile(j);
1325           ithSpectra=hSpectra.find(aTile->GetCellID());
1326           // double adc = aTile->GetPedestal()+aTile->GetIntegratedADC();
1327           double adc = aTile->GetIntegratedADC();
1328           double tot = aTile->GetIntegratedTOT();
1329           hspectraADCvsCellID->Fill(aTile->GetCellID(),adc);
1330           hspectraADCPedvsCellID->Fill(aTile->GetCellID(),aTile->GetPedestal());
1331           hspectraTOTvsCellID->Fill(aTile->GetCellID(),adc);
1332           if(ithSpectra!=hSpectra.end()){
1333             ithSpectra->second.FillExt(tot,adc, 0., 0.);
1334             ithSpectra->second.FillWaveform(aTile->GetADCWaveform());
1335           } else {
1336             RootOutputHist->cd("IndividualCells");
1337             hSpectra[aTile->GetCellID()]=TileSpectra("ped",2,aTile->GetCellID(),calib.GetTileCalib(aTile->GetCellID()),event.GetROtype(),debug);
1338             hSpectra[aTile->GetCellID()].FillExt(tot,adc, 0., 0.);
1339             hSpectra[aTile->GetCellID()].FillWaveform(aTile->GetADCWaveform());
1340             RootOutput->cd();
1341           }
1342         }
1343       }
1344     // }
1345     RootOutput->cd();
1346     TdataOut->Fill();
1347   }
1348   //RootOutput->cd();
1349   TdataOut->Write();
1350   TsetupIn->CloneTree()->Write();
1351 
1352   std::map<int,RunInfo>::iterator it=ri.find(runNr);
1353   TString outputDirPlots = GetPlotOutputDir();
1354   Double_t textSizeRel = 0.035;
1355 
1356   if (CalcBadChannel > 0 || ExtPlot > 0) {
1357     gSystem->Exec("mkdir -p "+outputDirPlots);
1358     StyleSettingsBasics("pdf");
1359     SetPlotStyle();  
1360   }
1361   
1362   if (CalcBadChannel > 0 ){
1363     //***********************************************************************************************
1364     //***** Monitoring histos for calib parameters & fits results of 1st iteration ******************
1365     //***********************************************************************************************
1366     int maxChannelPerLayer        = (setup->GetNMaxColumn()+1)*(setup->GetNMaxRow()+1);
1367     // monitoring applied pedestals
1368     TH1D* hBCvsCellID      = new TH1D( "hBC_vsCellID","Bad Channel vs CellID ; cell ID; BC flag ",
1369                                               setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5);
1370     hBCvsCellID->SetDirectory(0);
1371     TH2D* hBCVsLayer   = new TH2D( "hBCVsLayer","Bad Channel Map; layer; brd channel; BC flag  ",
1372                                               setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
1373     hBCVsLayer->SetDirectory(0);
1374 
1375     int currCells = 0;
1376     if ( debug > 0 && CalcBadChannel == 2)
1377       std::cout << "============================== starting bad channel extraction" << std::endl;
1378     if ( debug > 0 && CalcBadChannel == 1)
1379       std::cout << "============================== setting bad channel according to external map" << std::endl;
1380 
1381     for(ithSpectra=hSpectra.begin(); ithSpectra!=hSpectra.end(); ++ithSpectra){
1382       if (currCells%20 == 0 && currCells > 0 && debug > 0)
1383         std::cout << "============================== cell " <<  currCells << " / " << hSpectra.size() << " cells" << std::endl;
1384       currCells++;
1385       short bc   = 3;
1386       if (CalcBadChannel == 2){
1387         bc        = ithSpectra->second.DetermineBadChannel();
1388       } else if (CalcBadChannel == 1 && bcmap.size() > 0 ) {
1389         bcmapIte  = bcmap.find(ithSpectra->second.GetCellID());
1390         if(bcmapIte!=bcmap.end())
1391           bc        = bcmapIte->second;
1392         else 
1393           bc        = 3;
1394       } 
1395       
1396       long cellID   = ithSpectra->second.GetCellID();
1397       if (CalcBadChannel == 1)
1398         ithSpectra->second.SetBadChannelInCalib(bc);
1399       
1400       // initializing pedestal fits from calib file
1401       ithSpectra->second.InitializeNoiseFitsFromCalib();
1402       
1403       int layer     = setup->GetLayer(cellID);
1404       int chInLayer = setup->GetChannelInLayer(cellID);
1405       if (debug > 0 && bc > -1 && bc < 3)
1406         std::cout << "\t" << cellID << "\t" << layer << "\t" << setup->GetRow(cellID) << "\t" << setup->GetColumn(cellID)<< "\t" << setup->GetModule(cellID) << " - quality flag: " << bc << "\t" << calib.GetBadChannel(cellID) << "\t ped H: " << calib.GetPedestalMeanH(cellID) << "\t ped L: " << calib.GetPedestalMeanL(cellID)<< std::endl;
1407 
1408       hBCvsCellID->SetBinContent(hBCvsCellID->FindBin(cellID), calib.GetBadChannel(cellID));
1409       int bin2D     = hBCVsLayer->FindBin(layer,chInLayer);
1410       hBCVsLayer->SetBinContent(bin2D, calib.GetBadChannel(cellID));
1411     }
1412     hBCvsCellID->Write();
1413     hBCVsLayer->Write();
1414 
1415     //**********************************************************************
1416     // Create canvases for channel overview plotting
1417     //**********************************************************************
1418     // Get run info object
1419     // create directory for plot output
1420   
1421     TCanvas* canvas2DCorr = new TCanvas("canvasCorrPlots","",0,0,1450,1300);  // gives the page size
1422     DefaultCancasSettings( canvas2DCorr, 0.08, 0.13, 0.045, 0.07);
1423 
1424     canvas2DCorr->SetLogz(0);
1425     
1426     PlotSimple2DZRange( canvas2DCorr, hBCVsLayer, -10000, -10000, -0.1, 3.1, textSizeRel, Form("%s/BadChannelMap.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, "colz", true);    
1427     calib.SetBCCalib(true);
1428   }
1429     
1430   if (ExtPlot > 0){
1431     //***********************************************************************************************************
1432     //********************************* 8 Panel overview plot  **************************************************
1433     //***********************************************************************************************************
1434     //*****************************************************************
1435       // Test beam geometry (beam coming from viewer)
1436       //===========================================================
1437       //||    8 (4)    ||    7 (5)   ||    6 (6)   ||    5 (7)   ||  row 0
1438       //===========================================================
1439       //||    1 (0)    ||    2 (1)   ||    3 (2)   ||    4 (3)   ||  row 1
1440       //===========================================================
1441       //    col 0     col 1       col 2     col  3
1442       // rebuild pad geom in similar way (numbering -1)
1443     //*****************************************************************
1444     TCanvas* canvas8Panel;
1445     TPad* pad8Panel[8];
1446     Double_t topRCornerX[8];
1447     Double_t topRCornerY[8];
1448     Int_t textSizePixel = 30;
1449     Double_t relSize8P[8];
1450     CreateCanvasAndPadsFor8PannelTBPlot(canvas8Panel, pad8Panel,  topRCornerX, topRCornerY, relSize8P, textSizePixel);
1451 
1452     TCanvas* canvas8PanelProf;
1453     TPad* pad8PanelProf[8];
1454     Double_t topRCornerXProf[8];
1455     Double_t topRCornerYProf[8];
1456     Double_t relSize8PProf[8];
1457     CreateCanvasAndPadsFor8PannelTBPlot(canvas8PanelProf, pad8PanelProf,  topRCornerXProf, topRCornerYProf, relSize8PProf, textSizePixel, 0.045, "Prof", false);
1458 
1459     calib.PrintGlobalInfo();  
1460     // Double_t maxHG = ReturnMipPlotRangeDepVov(calib.GetVov(),true);
1461     // Double_t maxLG = ReturnMipPlotRangeDepVov(calib.GetVov(),false);
1462     Double_t maxHG = 600;
1463     Double_t maxLG = 200;
1464     
1465     std::cout << "plotting single layers" << std::endl;
1466     for (Int_t l = 0; l < setup->GetNMaxLayer()+1; l++){    
1467       if (l%10 == 0 && l > 0 && debug > 0)
1468         std::cout << "============================== layer " <<  l << " / " << setup->GetNMaxLayer() << " layers" << std::endl;     
1469       if (typeRO == ReadOut::Type::Caen) {
1470         PlotCorr2DFullLayer(canvas8PanelProf,pad8PanelProf, topRCornerXProf, topRCornerYProf, relSize8PProf,
1471                             textSizePixel, hSpectra, -20, 340, 4000, l, 0,
1472                             Form("%s/LGHG2D_Corr_Layer%02d.%s" ,outputDirPlots.Data(), l, plotSuffix.Data()), it->second);
1473       } else {
1474         PlotCorr2DFullLayer(canvas8PanelProf,pad8PanelProf, topRCornerXProf, topRCornerYProf, relSize8PProf,
1475                             textSizePixel, hSpectra, 0, it->second.samples+1, 1000, l, 0,
1476                             Form("%s/Waveform_Layer%02d.%s" ,outputDirPlots.Data(), l, plotSuffix.Data()), it->second);
1477         
1478       }
1479       if (ExtPlot > 1){
1480         PlotNoiseWithFitsFullLayer (canvas8Panel,pad8Panel, topRCornerX, topRCornerY, relSize8P, textSizePixel, 
1481                                     hSpectra, setup, true, -100, maxHG, 1.2, l, 0,
1482                                     Form("%s/SpectraWithNoiseFit_HG_Layer%02d.%s" ,outputDirPlots.Data(), l, plotSuffix.Data()), it->second);
1483         if (typeRO == ReadOut::Type::Caen){
1484           PlotNoiseWithFitsFullLayer (canvas8Panel,pad8Panel, topRCornerX, topRCornerY, relSize8P, textSizePixel, 
1485                                       hSpectra, setup, false, -20, maxLG, 1.2, l, 0,
1486                                       Form("%s/SpectraWithNoiseFit_LG_Layer%02d.%s" ,outputDirPlots.Data(), l, plotSuffix.Data()), it->second);
1487         }
1488       }
1489     }
1490     std::cout << "ending plotting single layers" << std::endl;
1491   }
1492   
1493   RootOutput->cd();
1494   
1495   std::cout<<"What is the value? <ped mean high>: "<<calib.GetAveragePedestalMeanHigh() << "\t good channels: " << calib.GetNumberOfChannelsWithBCflag(3) <<std::endl;
1496   if (IsCalibSaveToFile()){
1497     TString fileCalibPrint = RootOutputName;
1498     fileCalibPrint         = fileCalibPrint.ReplaceAll(".root","_calib.txt");
1499     calib.PrintCalibToFile(fileCalibPrint);
1500   }
1501   
1502   TcalibOut->Fill();
1503   TcalibOut->Write();
1504   
1505   RootOutput->Close();
1506   RootOutputHist->cd();
1507   if (typeRO == ReadOut::Type::Hgcroc){
1508      hspectraADCvsCellID->Write();
1509      hspectraADCPedvsCellID->Write();
1510      hspectraTOTvsCellID->Write();
1511   }
1512   RootOutputHist->cd("IndividualCells");
1513   for(ithSpectra=hSpectra.begin(); ithSpectra!=hSpectra.end(); ++ithSpectra){
1514     ithSpectra->second.WriteExt(true);
1515   }
1516 
1517 
1518   
1519   RootInput->Close();
1520   return true;
1521 }
1522 
1523 
1524 //***********************************************************************************************
1525 //*********************** intial scaling calculation function *********************************
1526 //***********************************************************************************************
1527 bool Analyses::GetScaling(void){
1528   std::cout<<"GetScaling"<<std::endl;
1529   
1530   std::map<int,RunInfo> ri=readRunInfosFromFile(RunListInputName.Data(),debug,0);
1531   
1532   std::map<int,TileSpectra> hSpectra;
1533   std::map<int,TileSpectra> hSpectraTrigg;
1534   std::map<int, TileSpectra>::iterator ithSpectra;
1535   std::map<int, TileSpectra>::iterator ithSpectraTrigg;
1536   
1537   ROOT::Math::MinimizerOptions::SetDefaultMinimizer("Minuit2", "Migrad");  
1538   
1539   std::cout << "Additional Output with histos being created: " << RootOutputNameHist.Data() << std::endl;
1540   if(Overwrite){
1541     std::cout << "recreating file with hists" << std::endl;
1542     RootOutputHist = new TFile(RootOutputNameHist.Data(),"RECREATE");
1543   } else{
1544     std::cout << "newly creating file with hists" << std::endl;
1545     RootOutputHist = new TFile(RootOutputNameHist.Data(),"CREATE");
1546   }
1547     
1548   //***********************************************************************************************
1549   //************************* first pass over tree to extract spectra *****************************
1550   //***********************************************************************************************
1551   // entering histoOutput file
1552   RootOutputHist->mkdir("IndividualCells");
1553   RootOutputHist->cd("IndividualCells");
1554 
1555   TcalibIn->GetEntry(0);
1556   // check whether calib should be overwritten based on external text file
1557   if (OverWriteCalib){
1558     calib.ReadCalibFromTextFile(ExternalCalibFile,debug);
1559   }
1560   
1561   int evts=TdataIn->GetEntries();
1562   int runNr = -1;
1563   if (maxEvents == -1){
1564     maxEvents = evts;
1565   } else {
1566     std::cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl;
1567     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;
1568     std::cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl;
1569   }
1570   ReadOut::Type typeRO = ReadOut::Type::Caen;
1571   int evtDeb = 5000;
1572   for(int i=0; i<evts && i < maxEvents; i++){
1573     TdataIn->GetEntry(i);
1574     if (i == 0){
1575       runNr = event.GetRunNumber();
1576       typeRO = event.GetROtype();
1577       if (typeRO != ReadOut::Type::Caen)
1578         evtDeb = 100;
1579       std::cout<< "Total number of events: " << evts << std::endl;
1580       std::cout<< "original run numbers calib: "<<calib.GetRunNumber() << "\t" << calib.GetRunNumberPed() << "\t" << calib.GetRunNumberMip() << std::endl;
1581       calib.SetRunNumberMip(runNr);
1582       calib.SetBeginRunTimeMip(event.GetBeginRunTimeAlt());
1583       std::cout<< "reset run numbers calib: "<< calib.GetRunNumber() << "\t" << calib.GetRunNumberPed() << "\t" << calib.GetRunNumberMip() << std::endl;
1584     }
1585     if (i%evtDeb == 0 && i > 0 && debug > 0) std::cout << "Reading " <<  i << " / " << evts << " events" << std::endl;
1586     if (i == 0 && debug > 2) std::cout << "decoding cell IDs" << std::endl ;
1587     for(int j=0; j<event.GetNTiles(); j++){
1588       
1589       if (typeRO == ReadOut::Type::Caen) {
1590         Caen* aTile=(Caen*)event.GetTile(j);
1591         if (i == 0 && debug > 2) std::cout << ((TString)setup->DecodeCellID(aTile->GetCellID())).Data() << std::endl;
1592         
1593         ithSpectra=hSpectra.find(aTile->GetCellID());
1594         double hgCorr = aTile->GetADCHigh()-calib.GetPedestalMeanH(aTile->GetCellID());
1595         double lgCorr = aTile->GetADCLow()-calib.GetPedestalMeanL(aTile->GetCellID());
1596             
1597         if(ithSpectra!=hSpectra.end()){
1598           ithSpectra->second.FillSpectra(lgCorr,hgCorr);
1599           if (hgCorr > 3*calib.GetPedestalSigH(aTile->GetCellID()) && lgCorr > 3*calib.GetPedestalSigL(aTile->GetCellID()) && hgCorr < 3900 )
1600             ithSpectra->second.FillCorr(lgCorr,hgCorr);
1601         } else {
1602           RootOutputHist->cd("IndividualCells");
1603           hSpectra[aTile->GetCellID()]=TileSpectra("mip1st",aTile->GetCellID(),calib.GetTileCalib(aTile->GetCellID()),event.GetROtype(), debug);
1604           hSpectra[aTile->GetCellID()].FillSpectra(lgCorr,hgCorr);;
1605           if (hgCorr > 3*calib.GetPedestalSigH(aTile->GetCellID()) && lgCorr > 3*calib.GetPedestalSigL(aTile->GetCellID() && hgCorr < 3900) )
1606             hSpectra[aTile->GetCellID()].FillCorr(lgCorr,hgCorr);;
1607           RootOutput->cd();
1608         }
1609       } else if (typeRO == ReadOut::Type::Hgcroc) {
1610         Hgcroc* aTile=(Hgcroc*)event.GetTile(j);
1611         if (i == 0 && debug > 2) std::cout << ((TString)setup->DecodeCellID(aTile->GetCellID())).Data() << std::endl;
1612         
1613         ithSpectra=hSpectra.find(aTile->GetCellID());
1614         double adc = aTile->GetIntegratedADC();
1615         double tot = aTile->GetIntegratedTOT();
1616             
1617         if(ithSpectra!=hSpectra.end()){
1618           ithSpectra->second.FillSpectra(tot,adc);
1619           ithSpectra->second.FillCorr(tot,adc);
1620         } else {
1621           RootOutputHist->cd("IndividualCells");
1622           hSpectra[aTile->GetCellID()]=TileSpectra("mip1st",aTile->GetCellID(),calib.GetTileCalib(aTile->GetCellID()),event.GetROtype(), debug);
1623           hSpectra[aTile->GetCellID()].FillSpectra(tot,adc);;
1624           hSpectra[aTile->GetCellID()].FillCorr(tot,adc);;
1625           RootOutput->cd();
1626         }
1627         
1628       }
1629     } 
1630     RootOutput->cd();
1631   }
1632   // TdataOut->Write();
1633   TsetupIn->CloneTree()->Write();
1634   
1635   RootOutputHist->cd();
1636  
1637   //***********************************************************************************************
1638   //***** Monitoring histos for calib parameters & fits results of 1st iteration ******************
1639   //***********************************************************************************************
1640   int maxChannelPerLayer        = (setup->GetNMaxColumn()+1)*(setup->GetNMaxRow()+1);
1641   // monitoring applied pedestals
1642   TH1D* hMeanPedHGvsCellID      = new TH1D( "hMeanPedHG_vsCellID","mean Ped High Gain vs CellID ; cell ID; #mu_{noise, HG} (arb. units) ",
1643                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5);
1644   hMeanPedHGvsCellID->SetDirectory(0);
1645   TH1D* hMeanPedHG              = new TH1D( "hMeanPedHG","mean Ped High Gain ; #mu_{noise, HG} (arb. units); counts ",
1646                                             500, -0.5, 500-0.5);
1647   hMeanPedHG->SetDirectory(0);
1648   TH2D* hspectraHGMeanVsLayer   = new TH2D( "hspectraHGMeanVsLayer","Mean Ped High Gain; layer; brd channel; #mu_{Ped HG} (arb. units) ",
1649                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
1650   hspectraHGMeanVsLayer->SetDirectory(0);
1651   TH2D* hspectraHGSigmaVsLayer  = new TH2D( "hspectraHGSigmaVsLayer","Sigma Ped High Gain; layer; brd channel; #sigma_{Ped HG} (arb. units) ",
1652                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
1653   hspectraHGSigmaVsLayer->SetDirectory(0);
1654   TH1D* hMeanPedLGvsCellID      = new TH1D( "hMeanPedLG_vsCellID","mean Ped Low Gain vs CellID ; cell ID; #mu_{noise, LG} (arb. units) ",
1655                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5);
1656   hMeanPedLGvsCellID->SetDirectory(0);
1657   TH1D* hMeanPedLG             = new TH1D( "hMeanPedLG","mean Ped Low Gain ; #mu_{noise, LG} (arb. units); counts ",
1658                                             500, -0.5, 500-0.5);
1659   hMeanPedLG->SetDirectory(0);
1660   TH2D* hspectraLGMeanVsLayer   = new TH2D( "hspectraLGMeanVsLayer","Mean Ped Low Gain; layer; brd channel; #mu_{PED LG} (arb. units) ",
1661                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
1662   hspectraLGMeanVsLayer->SetDirectory(0);
1663   TH2D* hspectraLGSigmaVsLayer  = new TH2D( "hspectraLGSigmaVsLayer","Sigma Ped Low Gain; layer; brd channel; #sigma_{Ped LG} (arb. units)",
1664                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
1665   hspectraLGSigmaVsLayer->SetDirectory(0);
1666   // monitoring 1st iteration mip fits
1667   TH2D* hspectraHGMaxVsLayer1st   = new TH2D( "hspectraHGMaxVsLayer_1st","Max High Gain; layer; brd channel; Max_{HG, 1^{st}} (arb. units) ",
1668                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
1669   hspectraHGMaxVsLayer1st->SetDirectory(0);
1670   TH2D* hspectraHGFWHMVsLayer1st   = new TH2D( "hspectraHGFWHMVsLayer_1st","FWHM High Gain; layer; brd channel; FWHM_{HG, 1^{st}} (arb. units) ",
1671                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
1672   hspectraHGFWHMVsLayer1st->SetDirectory(0);
1673   TH2D* hspectraHGLMPVVsLayer1st   = new TH2D( "hspectraHGLMPVVsLayer_1st","MPV High Gain; layer; brd channel; MPV_{HG, 1^{st}} (arb. units) ",
1674                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
1675   hspectraHGLMPVVsLayer1st->SetDirectory(0);
1676   TH2D* hspectraHGLSigmaVsLayer1st = new TH2D( "hspectraHGLSigmaVsLayer_1st","Sigma Landau High Gain; layer; brd channel; #sigma_{L,HG, 1^{st}} (arb. units) ",
1677                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
1678   hspectraHGLSigmaVsLayer1st->SetDirectory(0);
1679   TH2D* hspectraHGGSigmaVsLayer1st = new TH2D( "hspectraHGGSigmaVsLayer_1st","Sigma Gauss High Gain; layer; brd channel; #sigma_{G,HG, 1^{st}} (arb. units) ",
1680                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
1681   hspectraHGGSigmaVsLayer1st->SetDirectory(0);
1682   TH2D* hLGHGCorrVsLayer = new TH2D( "hLGHGCorrVsLayer","LG-HG corr; layer; brd channel; a_{LG-HG} (arb. units) ",
1683                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
1684   hLGHGCorrVsLayer->SetDirectory(0);
1685   TH2D* hHGLGCorrVsLayer = new TH2D( "hHGLGCorrVsLayer","HG-LG corr; layer; brd channel; a_{HG-LG} (arb. units) ",
1686                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
1687   hHGLGCorrVsLayer->SetDirectory(0);
1688   TH2D* hLGHGCorrOffsetVsLayer = new TH2D( "hLGHGCorrOffsetVsLayer","LG-HG corr offset; layer; brd channel; b_{LG-HG} (arb. units) ",
1689                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
1690   hLGHGCorrOffsetVsLayer->SetDirectory(0);
1691   TH2D* hHGLGCorrOffsetVsLayer = new TH2D( "hHGLGCorrVsLayer","HG-LG corr offset; layer; brd channel; b_{HG-LG} (arb. units) ",
1692                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
1693   hHGLGCorrOffsetVsLayer->SetDirectory(0);
1694   
1695   TH1D* hMaxHG1st             = new TH1D( "hMaxHG1st","Max High Gain ;Max_{HG, 1^{st}} (arb. units) ; counts ",
1696                                             2000, -0.5, 2000-0.5);
1697   hMaxHG1st->SetDirectory(0);
1698   TH1D* hLGHGCorr             = new TH1D( "hLGHGCorr","LG-HG corr ; a_{LG-HG} (arb. units) ; counts ",
1699                                             400, -20, 20);
1700   hLGHGCorr->SetDirectory(0);
1701   TH1D* hHGLGCorr             = new TH1D( "hHGLGCorr","LG-HG corr ; a_{HG-LG} (arb. units) ; counts ",
1702                                             400, -1., 1.);
1703   hHGLGCorr->SetDirectory(0);
1704 
1705   
1706   // fit pedestal
1707   double* parameters    = new double[6];
1708   double* parErrAndRes  = new double[6];
1709   bool isGood;
1710   int currCells = 0;
1711   if ( debug > 0)
1712     std::cout << "============================== starting fitting 1st iteration" << std::endl;
1713   
1714   for(ithSpectra=hSpectra.begin(); ithSpectra!=hSpectra.end(); ++ithSpectra){
1715     if (currCells%20 == 0 && currCells > 0 && debug > 0)
1716       std::cout << "============================== cell " <<  currCells << " / " << hSpectra.size() << " cells" << std::endl;
1717     currCells++;
1718     for (int p = 0; p < 6; p++){
1719       parameters[p]   = 0;
1720       parErrAndRes[p] = 0;
1721     }
1722     isGood        = ithSpectra->second.FitMipHG(parameters, parErrAndRes, debug, yearData, false, calib.GetVov(), 1);
1723     long cellID   = ithSpectra->second.GetCellID();
1724     int layer     = setup->GetLayer(cellID);
1725     int chInLayer = setup->GetChannelInLayer(cellID);
1726     
1727     // fill cross-check histos
1728     hMeanPedHGvsCellID->SetBinContent(hMeanPedHGvsCellID->FindBin(cellID), calib.GetPedestalMeanH(cellID));
1729     hMeanPedLGvsCellID->SetBinContent(hMeanPedLGvsCellID->FindBin(cellID), calib.GetPedestalMeanL(cellID));
1730     hMeanPedHG->Fill(calib.GetPedestalMeanH(cellID));
1731     hMeanPedLG->Fill(calib.GetPedestalMeanL(cellID));
1732     
1733     int bin2D     = hspectraHGMeanVsLayer->FindBin(layer,chInLayer);
1734     hspectraHGMeanVsLayer->SetBinContent(bin2D, calib.GetPedestalMeanH(cellID));
1735     hspectraHGSigmaVsLayer->SetBinContent(bin2D, calib.GetPedestalSigH(cellID));
1736     hspectraLGMeanVsLayer->SetBinContent(bin2D, calib.GetPedestalMeanL(cellID));
1737     hspectraLGSigmaVsLayer->SetBinContent(bin2D, calib.GetPedestalSigL(cellID));
1738 
1739     if (isGood){
1740       hspectraHGMaxVsLayer1st->SetBinContent(bin2D, parameters[4]);
1741       hspectraHGFWHMVsLayer1st->SetBinContent(bin2D, parameters[5]);
1742       hspectraHGLMPVVsLayer1st->SetBinContent(bin2D, parameters[1]);
1743       hspectraHGLMPVVsLayer1st->SetBinError(bin2D, parErrAndRes[1]);
1744       hspectraHGLSigmaVsLayer1st->SetBinContent(bin2D, parameters[0]);
1745       hspectraHGLSigmaVsLayer1st->SetBinError(bin2D, parErrAndRes[0]);
1746       hspectraHGGSigmaVsLayer1st->SetBinContent(bin2D, parameters[3]);
1747       hspectraHGGSigmaVsLayer1st->SetBinError(bin2D, parErrAndRes[3]);
1748       
1749       hMaxHG1st->Fill(parameters[4]);
1750     }
1751     
1752     isGood=ithSpectra->second.FitCorr(debug);
1753     if (ithSpectra->second.GetCorrModel(0)){
1754       hLGHGCorrVsLayer->SetBinContent(bin2D,ithSpectra->second.GetCorrModel(0)->GetParameter(1));
1755       hLGHGCorrVsLayer->SetBinError(bin2D,ithSpectra->second.GetCorrModel(0)->GetParError(1));
1756       hLGHGCorrOffsetVsLayer->SetBinContent(bin2D,ithSpectra->second.GetCorrModel(0)->GetParameter(0));
1757       hLGHGCorrOffsetVsLayer->SetBinError(bin2D,ithSpectra->second.GetCorrModel(0)->GetParError(0));
1758       hLGHGCorr->Fill(ithSpectra->second.GetCorrModel(0)->GetParameter(1));
1759     } 
1760     if (ithSpectra->second.GetCorrModel(1)){
1761       hHGLGCorrVsLayer->SetBinContent(bin2D,ithSpectra->second.GetCorrModel(1)->GetParameter(1));
1762       hHGLGCorrVsLayer->SetBinError(bin2D,ithSpectra->second.GetCorrModel(1)->GetParError(1));    
1763       hHGLGCorrOffsetVsLayer->SetBinContent(bin2D,ithSpectra->second.GetCorrModel(1)->GetParameter(0));
1764       hHGLGCorrOffsetVsLayer->SetBinError(bin2D,ithSpectra->second.GetCorrModel(1)->GetParError(0));    
1765       hHGLGCorr->Fill(ithSpectra->second.GetCorrModel(1)->GetParameter(1));
1766     }
1767   }
1768   if ( debug > 0)
1769     std::cout << "============================== done fitting 1st iteration" << std::endl;
1770 
1771   TH1D* hHGTileSum[20];
1772   for (int c = 0; c < maxChannelPerLayer; c++ ){
1773     hHGTileSum[c] = new TH1D( Form("hHGTileSum_absChB%d",c),"av ADC surrounding cells ; ADC (arb. units); counts ",
1774                               4000, -0.5, 4000-0.5);
1775     hHGTileSum[c]->SetDirectory(0);
1776   }
1777   
1778   // setup trigger sel
1779   TRandom3* rand    = new TRandom3();
1780   Int_t localTriggerTiles = 4;
1781   double factorMinTrigg   = 0.5;
1782   double factorMaxTrigg   = 4.;
1783   if (yearData == 2023){
1784     localTriggerTiles = 6;
1785     factorMaxTrigg    = 3.;
1786   }
1787   RootOutputHist->mkdir("IndividualCellsTrigg");
1788   RootOutputHist->cd("IndividualCellsTrigg");  
1789   //***********************************************************************************************
1790   //************************* first pass over tree to extract spectra *****************************
1791   //***********************************************************************************************  
1792   int actCh1st = 0;
1793   double averageScale = calib.GetAverageScaleHigh(actCh1st);
1794   double meanFWHMHG   = calib.GetAverageScaleWidthHigh();
1795   double avHGLGCorr   = calib.GetAverageHGLGCorr();
1796   double avHGLGOffCorr= calib.GetAverageHGLGCorrOff();
1797   double avLGHGCorr   = calib.GetAverageLGHGCorr();
1798   double avLGHGOffCorr= calib.GetAverageLGHGCorrOff();
1799   std::cout << "average HG mip: " << averageScale << "\t active ch: "<< actCh1st<< std::endl;
1800   for(int i=0; i<evts && i < maxEvents; i++){
1801     TdataIn->GetEntry(i);
1802     if (i%evtDeb == 0 && i > 0 && debug > 0) std::cout << "Reading " <<  i << " / " << evts << " events" << std::endl;
1803     for(int j=0; j<event.GetNTiles(); j++){
1804       if (typeRO == ReadOut::Type::Caen) {
1805         Caen* aTile=(Caen*)event.GetTile(j);
1806         if (i == 0 && debug > 2) std::cout << ((TString)setup->DecodeCellID(aTile->GetCellID())).Data() << std::endl;
1807         long currCellID = aTile->GetCellID();
1808         
1809         // read current tile
1810         ithSpectraTrigg=hSpectraTrigg.find(aTile->GetCellID());
1811         double hgCorr = aTile->GetADCHigh()-calib.GetPedestalMeanH(aTile->GetCellID());
1812         double lgCorr = aTile->GetADCLow()-calib.GetPedestalMeanL(aTile->GetCellID());
1813 
1814         aTile->SetLocalTriggerPrimitive(event.CalculateLocalMuonTrigg(calib, rand, currCellID, localTriggerTiles, avLGHGCorr));
1815         // estimate local muon trigger
1816         bool localMuonTrigg = event.InspectIfLocalMuonTrigg(currCellID, averageScale, factorMinTrigg, factorMaxTrigg);
1817         int chInLayer = setup->GetChannelInLayer(currCellID);    
1818         hHGTileSum[chInLayer]->Fill(aTile->GetLocalTriggerPrimitive());
1819         
1820         if(ithSpectraTrigg!=hSpectraTrigg.end()){
1821           ithSpectraTrigg->second.FillTrigger(aTile->GetLocalTriggerPrimitive());
1822         } else {
1823           RootOutputHist->cd("IndividualCellsTrigg");
1824           hSpectraTrigg[currCellID]=TileSpectra("mipTrigg",currCellID,calib.GetTileCalib(currCellID),event.GetROtype(),debug);
1825           hSpectraTrigg[currCellID].FillTrigger(aTile->GetLocalTriggerPrimitive());;
1826           RootOutput->cd();
1827         }
1828         // only fill tile spectra if 4 surrounding tiles on average are compatible with muon response
1829         if (localMuonTrigg){
1830           aTile->SetLocalTriggerBit(1);
1831           ithSpectraTrigg=hSpectraTrigg.find(aTile->GetCellID());
1832           ithSpectraTrigg->second.FillSpectra(lgCorr,hgCorr);
1833           if (hgCorr > 3*calib.GetPedestalSigH(currCellID) && lgCorr > 3*calib.GetPedestalSigL(currCellID) && hgCorr < 3900 )
1834             ithSpectraTrigg->second.FillCorr(lgCorr,hgCorr);
1835         }
1836       } else if (typeRO == ReadOut::Type::Hgcroc) {
1837         Hgcroc* aTile=(Hgcroc*)event.GetTile(j);
1838         if (i == 0 && debug > 2) std::cout << ((TString)setup->DecodeCellID(aTile->GetCellID())).Data() << std::endl;
1839         long currCellID = aTile->GetCellID();
1840         
1841         ithSpectraTrigg=hSpectraTrigg.find(aTile->GetCellID());
1842         // double adc = aTile->GetPedestal()+aTile->GetIntegratedADC();
1843         double adc = aTile->GetIntegratedADC();
1844         double tot = aTile->GetIntegratedTOT();
1845             
1846         aTile->SetLocalTriggerPrimitive(event.CalculateLocalMuonTrigg(calib, rand, currCellID, localTriggerTiles, 0));
1847         // estimate local muon trigger
1848         bool localMuonTrigg = event.InspectIfLocalMuonTrigg(currCellID, 0, factorMinTrigg, factorMaxTrigg);
1849         int chInLayer = setup->GetChannelInLayer(currCellID);    
1850         hHGTileSum[chInLayer]->Fill(aTile->GetLocalTriggerPrimitive());
1851         
1852         if(ithSpectraTrigg!=hSpectraTrigg.end()){
1853           ithSpectraTrigg->second.FillTrigger(aTile->GetLocalTriggerPrimitive());
1854         } else {
1855           RootOutputHist->cd("IndividualCellsTrigg");
1856           hSpectraTrigg[currCellID]=TileSpectra("mipTrigg",currCellID,calib.GetTileCalib(currCellID),event.GetROtype(),debug);
1857           hSpectraTrigg[currCellID].FillTrigger(aTile->GetLocalTriggerPrimitive());;
1858           RootOutput->cd();
1859         }
1860         // only fill tile spectra if 4 surrounding tiles on average are compatible with muon response
1861         if (localMuonTrigg){
1862           aTile->SetLocalTriggerBit(1);
1863           ithSpectraTrigg=hSpectraTrigg.find(aTile->GetCellID());
1864           ithSpectraTrigg->second.FillSpectra(tot,adc);
1865           ithSpectraTrigg->second.FillCorr(tot,adc);
1866         }
1867       }
1868     }
1869     RootOutput->cd();
1870     TdataOut->Fill();
1871   }
1872   TdataOut->Write();
1873   
1874   //***********************************************************************************************
1875   //***** Monitoring histos for fits results of 2nd iteration ******************
1876   //***********************************************************************************************
1877   RootOutputHist->cd();
1878   
1879   // monitoring trigger 
1880   TH2D* hmipTriggers              = new TH2D( "hmipTriggers","muon triggers; layer; brd channel; counts ",
1881                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
1882   hmipTriggers->SetDirectory(0);
1883   TH2D* hSuppresionNoise          = new TH2D( "hSuppresionNoise","S/B noise region; layer; brd channel; S/B noise region",
1884                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
1885   hSuppresionNoise->SetDirectory(0);
1886   TH2D* hSuppresionSignal         = new TH2D( "hSuppresionSignal","S/B signal region; layer; brd channel; S/B signal region",
1887                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
1888   hSuppresionSignal->SetDirectory(0);
1889 
1890   // monitoring 2nd iteration mip fits
1891   TH2D* hspectraHGMaxVsLayer2nd   = new TH2D( "hspectraHGMaxVsLayer_2nd","Max High Gain; layer; brd channel; Max_{HG, 2^{nd}} (arb. units) ",
1892                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
1893   hspectraHGMaxVsLayer2nd->SetDirectory(0);
1894   TH2D* hspectraHGFWHMVsLayer2nd   = new TH2D( "hspectraHGFWHMVsLayer_2nd","FWHM High Gain; layer; brd channel; FWHM_{HG, 2^{nd}} (arb. units) ",
1895                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
1896   hspectraHGFWHMVsLayer2nd->SetDirectory(0);
1897   TH2D* hspectraHGLMPVVsLayer2nd   = new TH2D( "hspectraHGLMPVVsLayer_2nd","MPV High Gain; layer; brd channel; MPV_{HG, 2^{nd}} (arb. units) ",
1898                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
1899   hspectraHGLMPVVsLayer2nd->SetDirectory(0);
1900   TH2D* hspectraHGLSigmaVsLayer2nd = new TH2D( "hspectraHGLSigmaVsLayer_2nd","Sigma Landau High Gain; layer; brd channel; #sigma_{L,HG, 2^{nd}} (arb. units) ",
1901                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
1902   hspectraHGLSigmaVsLayer2nd->SetDirectory(0);
1903   TH2D* hspectraHGGSigmaVsLayer2nd = new TH2D( "hspectraHGGSigmaVsLayer_2nd","Sigma Gauss High Gain; layer; brd channel; #sigma_{G,HG, 2^{nd}} (arb. units) ",
1904                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
1905   hspectraHGGSigmaVsLayer2nd->SetDirectory(0);
1906   TH2D* hspectraLGMaxVsLayer2nd   = new TH2D( "hspectraLGMaxVsLayer_2nd","Max High Gain; layer; brd channel; Max_{LG, 2^{nd}} (arb. units) ",
1907                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
1908   hspectraLGMaxVsLayer2nd->SetDirectory(0);
1909   TH2D* hspectraLGFWHMVsLayer2nd   = new TH2D( "hspectraLGFWHMVsLayer_2nd","FWHM High Gain; layer; brd channel; FWHM_{LG, 2^{nd}} (arb. units) ",
1910                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
1911   hspectraLGFWHMVsLayer2nd->SetDirectory(0);
1912   TH2D* hspectraLGLMPVVsLayer2nd   = new TH2D( "hspectraLGLMPVVsLayer_2nd","MPV High Gain; layer; brd channel; MPV_{LG, 2^{nd}} (arb. units) ",
1913                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
1914   hspectraLGLMPVVsLayer2nd->SetDirectory(0);
1915   TH2D* hspectraLGLSigmaVsLayer2nd = new TH2D( "hspectraLGLSigmaVsLayer_2nd","Sigma Landau High Gain; layer; brd channel; #sigma_{L,LG, 2^{nd}} (arb. units) ",
1916                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
1917   hspectraLGLSigmaVsLayer2nd->SetDirectory(0);
1918   TH2D* hspectraLGGSigmaVsLayer2nd = new TH2D( "hspectraLGGSigmaVsLayer_2nd","Sigma Gauss High Gain; layer; brd channel; #sigma_{G,LG, 2^{nd}} (arb. units) ",
1919                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
1920   hspectraLGGSigmaVsLayer2nd->SetDirectory(0);
1921 
1922   TH1D* hMaxHG2nd             = new TH1D( "hMaxHG2nd","Max High Gain ;Max_{HG, 2^{nd}} (arb. units) ; counts ",
1923                                             2000, -0.5, 2000-0.5);
1924   hMaxHG2nd->SetDirectory(0);
1925   TH1D* hMaxLG2nd             = new TH1D( "hMaxLG2nd","Max High Gain ;Max_{LG, 2^{nd}} (arb. units) ; counts ",
1926                                             400, -0.5, 400-0.5);
1927   hMaxLG2nd->SetDirectory(0);
1928 
1929 
1930   currCells = 0;
1931   if ( debug > 0)
1932     std::cout << "============================== starting fitting 2nd iteration" << std::endl;
1933   for(ithSpectraTrigg=hSpectraTrigg.begin(); ithSpectraTrigg!=hSpectraTrigg.end(); ++ithSpectraTrigg){
1934     if (currCells%20 == 0 && currCells > 0 && debug > 0)
1935       std::cout << "============================== cell " <<  currCells << " / " << hSpectraTrigg.size() << " cells" << std::endl;
1936     currCells++;
1937     for (int p = 0; p < 6; p++){
1938       parameters[p]   = 0;
1939       parErrAndRes[p] = 0;
1940     }
1941     isGood=ithSpectraTrigg->second.FitMipHG(parameters, parErrAndRes, debug, yearData, false, calib.GetVov(), averageScale);
1942     
1943     long cellID     = ithSpectraTrigg->second.GetCellID();
1944     int layer       = setup->GetLayer(cellID);
1945     int chInLayer   = setup->GetChannelInLayer(cellID);    
1946     int bin2D       = hspectraHGMeanVsLayer->FindBin(layer,chInLayer);
1947 
1948     double pedSigHG = 0;
1949     double maxBin   = 0;
1950     if (typeRO == ReadOut::Type::Caen){
1951       pedSigHG = calib.GetPedestalSigH(cellID);
1952       maxBin   = 3800;
1953     } else {
1954       pedSigHG = 20;
1955       maxBin   = 1024;
1956     }
1957     Int_t binNLow   = ithSpectraTrigg->second.GetHG()->FindBin(-1*pedSigHG);
1958     Int_t binNHigh  = ithSpectraTrigg->second.GetHG()->FindBin(3*pedSigHG);
1959     Int_t binSHigh  = ithSpectraTrigg->second.GetHG()->FindBin(maxBin);
1960     
1961     double S_NoiseR = ithSpectraTrigg->second.GetHG()->Integral(binNLow, binNHigh);
1962     double S_SigR   = ithSpectraTrigg->second.GetHG()->Integral(binNHigh, binSHigh);
1963     
1964     ithSpectra      = hSpectra.find(cellID);
1965     double B_NoiseR = ithSpectra->second.GetHG()->Integral(binNLow , binNHigh);
1966     double B_SigR   = ithSpectra->second.GetHG()->Integral(binNHigh, binSHigh);
1967     
1968     double SB_NoiseR  = (B_NoiseR != 0.) ? S_NoiseR/B_NoiseR : 0;
1969     double SB_SigR    = (B_SigR != 0.) ? S_SigR/B_SigR : 0;
1970     
1971     hmipTriggers->SetBinContent(bin2D, ithSpectraTrigg->second.GetHG()->GetEntries());
1972     hSuppresionNoise->SetBinContent(bin2D, SB_NoiseR);
1973     hSuppresionSignal->SetBinContent(bin2D, SB_SigR);
1974     if (isGood){
1975       hspectraHGMaxVsLayer2nd->SetBinContent(bin2D, parameters[4]);
1976       hspectraHGFWHMVsLayer2nd->SetBinContent(bin2D, parameters[5]);
1977       hspectraHGLMPVVsLayer2nd->SetBinContent(bin2D, parameters[1]);
1978       hspectraHGLMPVVsLayer2nd->SetBinError(bin2D, parErrAndRes[1]);
1979       hspectraHGLSigmaVsLayer2nd->SetBinContent(bin2D, parameters[0]);
1980       hspectraHGLSigmaVsLayer2nd->SetBinError(bin2D, parErrAndRes[0]);
1981       hspectraHGGSigmaVsLayer2nd->SetBinContent(bin2D, parameters[3]);
1982       hspectraHGGSigmaVsLayer2nd->SetBinError(bin2D, parErrAndRes[3]);
1983       hMaxHG2nd->Fill(parameters[4]);
1984     }
1985     for (int p = 0; p < 6; p++){
1986       parameters[p]   = 0;
1987       parErrAndRes[p] = 0;
1988     }
1989     isGood=ithSpectraTrigg->second.FitMipLG(parameters, parErrAndRes, debug, yearData, false, 1);
1990     if (isGood){
1991       hspectraLGMaxVsLayer2nd->SetBinContent(bin2D, parameters[4]);
1992       hspectraLGFWHMVsLayer2nd->SetBinContent(bin2D, parameters[5]);
1993       hspectraLGLMPVVsLayer2nd->SetBinContent(bin2D, parameters[1]);
1994       hspectraLGLMPVVsLayer2nd->SetBinError(bin2D, parErrAndRes[1]);
1995       hspectraLGLSigmaVsLayer2nd->SetBinContent(bin2D, parameters[0]);
1996       hspectraLGLSigmaVsLayer2nd->SetBinError(bin2D, parErrAndRes[0]);
1997       hspectraLGGSigmaVsLayer2nd->SetBinContent(bin2D, parameters[3]);
1998       hspectraLGGSigmaVsLayer2nd->SetBinError(bin2D, parErrAndRes[3]);
1999       hMaxLG2nd->Fill(parameters[4]);
2000     }
2001   }
2002   if ( debug > 0)
2003     std::cout << "============================== done fitting 2nd iteration" << std::endl;
2004   int actCh2nd = 0;
2005   double averageScaleUpdated    = calib.GetAverageScaleHigh(actCh2nd);
2006   double meanFWHMHGUpdated      = calib.GetAverageScaleWidthHigh();
2007   double averageScaleLowUpdated = calib.GetAverageScaleLow();
2008   double meanFWHMLGUpdated      = calib.GetAverageScaleWidthLow();
2009   std::cout << "average 1st iteration HG mip: " << averageScale << "\t act. channels: " <<   actCh1st << std::endl;
2010   std::cout << "average 2nd iteration  HG mip: " << averageScaleUpdated << "\t LG mip: " << averageScaleLowUpdated << "\t act. channels: " <<   actCh2nd << std::endl;
2011   
2012   RootOutput->cd();
2013   
2014   if (IsCalibSaveToFile()){
2015     TString fileCalibPrint = RootOutputName;
2016     fileCalibPrint         = fileCalibPrint.ReplaceAll(".root","_calib.txt");
2017     calib.PrintCalibToFile(fileCalibPrint);
2018   }
2019 
2020   TcalibOut->Fill();
2021   TcalibOut->Write();
2022   RootOutput->Close();
2023 
2024 
2025   RootOutputHist->cd("IndividualCells");
2026     for(ithSpectra=hSpectra.begin(); ithSpectra!=hSpectra.end(); ++ithSpectra){
2027       ithSpectra->second.Write(true);
2028     }
2029   RootOutputHist->cd("IndividualCellsTrigg");
2030     for(ithSpectra=hSpectraTrigg.begin(); ithSpectra!=hSpectraTrigg.end(); ++ithSpectra){
2031       ithSpectra->second.Write(true);
2032     }
2033   RootOutputHist->cd();
2034     hMeanPedHGvsCellID->Write();
2035     hMeanPedHG->Write();
2036     hMeanPedLGvsCellID->Write();
2037     hMeanPedLG->Write();
2038     
2039     hspectraHGMeanVsLayer->Write();
2040     hspectraLGMeanVsLayer->Write();
2041     hspectraHGSigmaVsLayer->Write();
2042     hspectraLGSigmaVsLayer->Write();
2043     hspectraHGMaxVsLayer1st->Write();
2044     hspectraHGFWHMVsLayer1st->Write();
2045     hspectraHGLMPVVsLayer1st->Write();
2046     hspectraHGLSigmaVsLayer1st->Write();
2047     hspectraHGGSigmaVsLayer1st->Write();
2048     hMaxHG1st->Write();
2049     
2050     hLGHGCorrVsLayer->Write();
2051     hLGHGCorrOffsetVsLayer->Write();
2052     hHGLGCorrVsLayer->Write();
2053     hHGLGCorrOffsetVsLayer->Write();
2054     hLGHGCorr->Write();
2055     hHGLGCorr->Write();
2056     
2057     hspectraHGMaxVsLayer2nd->Write();
2058     hspectraHGFWHMVsLayer2nd->Write();
2059     hspectraHGLMPVVsLayer2nd->Write();
2060     hspectraHGLSigmaVsLayer2nd->Write();
2061     hspectraHGGSigmaVsLayer2nd->Write();
2062     hMaxHG2nd->Write();
2063     
2064     hspectraLGMaxVsLayer2nd->Write();
2065     hspectraLGFWHMVsLayer2nd->Write();
2066     hspectraLGLMPVVsLayer2nd->Write();
2067     hspectraLGLSigmaVsLayer2nd->Write();
2068     hspectraLGGSigmaVsLayer2nd->Write();
2069     hMaxLG2nd->Write();
2070     
2071     for (int c = 0; c< maxChannelPerLayer; c++)
2072       hHGTileSum[c]->Write();
2073     hmipTriggers->Write();
2074     hSuppresionNoise->Write();
2075     hSuppresionSignal->Write();
2076   // fill calib tree & write it
2077   // close open root files
2078   RootOutputHist->Write();
2079   RootOutputHist->Close();
2080 
2081   RootInput->Close();
2082 
2083   // Get run info object
2084   std::map<int,RunInfo>::iterator it=ri.find(runNr);
2085   
2086   // create directory for plot output
2087   TString outputDirPlots = GetPlotOutputDir();
2088   gSystem->Exec("mkdir -p "+outputDirPlots);
2089   
2090   //**********************************************************************
2091   // Create canvases for channel overview plotting
2092   //**********************************************************************
2093   Double_t textSizeRel = 0.035;
2094   StyleSettingsBasics("pdf");
2095   SetPlotStyle();
2096 
2097   TCanvas* canvas2DCorr = new TCanvas("canvasCorrPlots","",0,0,1450,1300);  // gives the page size
2098   DefaultCancasSettings( canvas2DCorr, 0.08, 0.13, 0.045, 0.07);
2099 
2100   canvas2DCorr->SetLogz(0);
2101   PlotSimple2D( canvas2DCorr, hspectraHGMeanVsLayer, -10000, -10000, textSizeRel, Form("%s/HG_NoiseMean.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
2102   PlotSimple2D( canvas2DCorr, hspectraHGSigmaVsLayer,-10000, -10000, textSizeRel, Form("%s/HG_NoiseSigma.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1,  kFALSE, "colz", true);
2103   PlotSimple2D( canvas2DCorr, hspectraHGMaxVsLayer1st, -10000, -10000, textSizeRel, Form("%s/HG_MaxMip_1st.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true, Form( "#LT Max_{HG} #GT = %.1f", averageScale));
2104   PlotSimple2D( canvas2DCorr, hspectraHGFWHMVsLayer1st, -10000, -10000, textSizeRel, Form("%s/HG_FWHMMip_1st.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true, Form( "#LT FWHM_{HG} #GT = %.1f", meanFWHMHG));
2105   PlotSimple2D( canvas2DCorr, hspectraHGLMPVVsLayer1st, -10000, -10000, textSizeRel, Form("%s/HG_LandMPVMip_1st.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
2106   PlotSimple2D( canvas2DCorr, hspectraHGLSigmaVsLayer1st, -10000, -10000, textSizeRel, Form("%s/HG_LandSigMip_1st.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
2107   PlotSimple2D( canvas2DCorr, hspectraHGGSigmaVsLayer1st, -10000, -10000, textSizeRel, Form("%s/HG_GaussSigMip_1st.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
2108   PlotSimple2D( canvas2DCorr, hspectraHGMaxVsLayer2nd, -10000, -10000, textSizeRel, Form("%s/HG_MaxMip_2nd.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true, Form( "#LT Max_{HG} #GT = %.1f", averageScaleUpdated));
2109   PlotSimple2D( canvas2DCorr, hspectraHGFWHMVsLayer2nd, -10000, -10000, textSizeRel, Form("%s/HG_FWHMMip_2nd.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true, Form( "#LT FWHM_{HG} #GT = %.1f", meanFWHMHGUpdated));
2110   PlotSimple2D( canvas2DCorr, hspectraHGLMPVVsLayer2nd, -10000, -10000, textSizeRel, Form("%s/HG_LandMPVMip_2nd.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
2111   PlotSimple2D( canvas2DCorr, hspectraHGLSigmaVsLayer2nd, -10000, -10000, textSizeRel, Form("%s/HG_LandSigMip_2nd.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
2112   PlotSimple2D( canvas2DCorr, hspectraHGGSigmaVsLayer2nd, -10000, -10000, textSizeRel, Form("%s/HG_GaussSigMip_2nd.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
2113   canvas2DCorr->SetLogz(1);
2114   TString drawOpt = "colz";
2115   if (yearData == 2023)
2116     drawOpt = "colz,text";
2117   PlotSimple2D( canvas2DCorr, hmipTriggers, -10000, -10000, textSizeRel, Form("%s/MuonTriggers.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, drawOpt, true, Form( "evt. collected = %d", evts));
2118   PlotSimple2D( canvas2DCorr, hSuppresionNoise, -10000, -10000, textSizeRel, Form("%s/SuppressionNoise.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, drawOpt, true);
2119   PlotSimple2D( canvas2DCorr, hSuppresionSignal, -10000, -10000, textSizeRel, Form("%s/SuppressionSignal.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, drawOpt, true);
2120   
2121   canvas2DCorr->SetLogz(0);
2122   PlotSimple2D( canvas2DCorr, hspectraLGMeanVsLayer, -10000, -10000, textSizeRel, Form("%s/LG_NoiseMean.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
2123   PlotSimple2D( canvas2DCorr, hspectraLGSigmaVsLayer, -10000, -10000, textSizeRel, Form("%s/LG_NoiseSigma.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
2124   PlotSimple2D( canvas2DCorr, hspectraLGMaxVsLayer2nd, -10000, -10000, textSizeRel, Form("%s/LG_MaxMip_2nd.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true, Form( "#LT Max_{LG} #GT = %.1f", averageScaleLowUpdated));
2125   PlotSimple2D( canvas2DCorr, hspectraLGFWHMVsLayer2nd, -10000, -10000, textSizeRel, Form("%s/LG_FWHMMip_2nd.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true, Form( "#LT FWHM_{LG} #GT = %.1f", meanFWHMLGUpdated));
2126   PlotSimple2D( canvas2DCorr, hspectraLGLMPVVsLayer2nd, -10000, -10000, textSizeRel, Form("%s/LG_LandMPVMip_2nd.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
2127   PlotSimple2D( canvas2DCorr, hspectraLGLSigmaVsLayer2nd, -10000, -10000, textSizeRel, Form("%s/LG_LandSigMip_2nd.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
2128   PlotSimple2D( canvas2DCorr, hspectraLGGSigmaVsLayer2nd, -10000, -10000, textSizeRel, Form("%s/LG_GaussSigMip_2nd.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
2129 
2130   PlotSimple2D( canvas2DCorr, hLGHGCorrVsLayer, -10000, -10000, textSizeRel, Form("%s/LG_HG_Corr.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kTRUE, "colz", true, Form( "#LT a_{LGHG} #GT = %.1f", avLGHGCorr));
2131   PlotSimple2D( canvas2DCorr, hLGHGCorrOffsetVsLayer, -10000, -10000, textSizeRel, Form("%s/LG_HG_CorrOffset.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kTRUE, "colz", true, Form( "#LT b_{LGHG} #GT = %.1f", avLGHGOffCorr));
2132   PlotSimple2D( canvas2DCorr, hHGLGCorrVsLayer, -10000, -10000, textSizeRel, Form("%s/HG_LG_Corr.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kTRUE, "colz", true, Form( "#LT a_{HGLG} #GT = %.1f", avHGLGCorr));
2133   PlotSimple2D( canvas2DCorr, hHGLGCorrOffsetVsLayer, -10000, -10000, textSizeRel, Form("%s/HG_LG_CorrOffset.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kTRUE, "colz", true, Form( "#LT b_{HGLG} #GT = %.1f", avHGLGOffCorr));
2134   
2135   if (ExtPlot > 0){
2136     //***********************************************************************************************************
2137     //********************************* 8 Panel overview plot  **************************************************
2138     //***********************************************************************************************************
2139     //*****************************************************************
2140       // Test beam geometry (beam coming from viewer)
2141       //===========================================================
2142       //||    8 (4)    ||    7 (5)   ||    6 (6)   ||    5 (7)   ||  row 0
2143       //===========================================================
2144       //||    1 (0)    ||    2 (1)   ||    3 (2)   ||    4 (3)   ||  row 1
2145       //===========================================================
2146       //    col 0     col 1       col 2     col  3
2147       // rebuild pad geom in similar way (numbering -1)
2148     //*****************************************************************
2149     TCanvas* canvas8Panel;
2150     TPad* pad8Panel[8];
2151     Double_t topRCornerX[8];
2152     Double_t topRCornerY[8];
2153     Int_t textSizePixel = 30;
2154     Double_t relSize8P[8];
2155     CreateCanvasAndPadsFor8PannelTBPlot(canvas8Panel, pad8Panel,  topRCornerX, topRCornerY, relSize8P, textSizePixel);
2156 
2157     TCanvas* canvas8PanelProf;
2158     TPad* pad8PanelProf[8];
2159     Double_t topRCornerXProf[8];
2160     Double_t topRCornerYProf[8];
2161     Double_t relSize8PProf[8];
2162     CreateCanvasAndPadsFor8PannelTBPlot(canvas8PanelProf, pad8PanelProf,  topRCornerXProf, topRCornerYProf, relSize8PProf, textSizePixel, 0.045, "Prof", false);
2163 
2164     calib.PrintGlobalInfo();  
2165     Double_t maxHG = ReturnMipPlotRangeDepVov(calib.GetVov(),true);
2166     Double_t maxLG = ReturnMipPlotRangeDepVov(calib.GetVov(),false);
2167     std::cout << "plotting single layers" << std::endl;
2168 
2169     for (Int_t l = 0; l < setup->GetNMaxLayer()+1; l++){    
2170       if (l%10 == 0 && l > 0 && debug > 0)
2171         std::cout << "============================== layer " <<  l << " / " << setup->GetNMaxLayer() << " layers" << std::endl;
2172       if (ExtPlot > 0){
2173         PlotMipWithFitsFullLayer (canvas8Panel,pad8Panel, topRCornerX, topRCornerY, relSize8P, textSizePixel, 
2174                                   hSpectra, hSpectraTrigg, setup, true, -100, maxHG, 1.2, l, 0,
2175                                   Form("%s/MIP_HG_Layer%02d.%s" ,outputDirPlots.Data(), l, plotSuffix.Data()), it->second);
2176         PlotTriggerPrimWithFitsFullLayer (canvas8Panel,pad8Panel, topRCornerX, topRCornerY, relSize8P, textSizePixel, 
2177                                           hSpectraTrigg, setup, averageScale, factorMinTrigg, factorMaxTrigg,
2178                                           0, maxHG*2, 1.2, l, 0, Form("%s/TriggPrimitive_Layer%02d.%s" ,outputDirPlots.Data(), l, plotSuffix.Data()), it->second);
2179         PlotCorrWithFitsFullLayer(canvas8PanelProf,pad8PanelProf, topRCornerXProf, topRCornerYProf, relSize8PProf, textSizePixel, 
2180                                   hSpectra, 0, -20, 800, 3900, l, 0,
2181                                   Form("%s/LGHG_Corr_Layer%02d.%s" ,outputDirPlots.Data(), l, plotSuffix.Data()), it->second);
2182       }
2183       if (ExtPlot > 1){
2184         PlotMipWithFitsFullLayer (canvas8Panel,pad8Panel, topRCornerX, topRCornerY, relSize8P, textSizePixel, 
2185                                   hSpectra, hSpectraTrigg, setup, false, -30, maxLG, 1.2, l, 0,
2186                                   Form("%s/MIP_LG_Layer%02d.%s" ,outputDirPlots.Data(), l, plotSuffix.Data()), it->second);
2187         PlotCorrWithFitsFullLayer(canvas8PanelProf,pad8PanelProf, topRCornerXProf, topRCornerYProf, relSize8PProf, textSizePixel, 
2188                                   hSpectra, 1, -100, 4000, 340, l, 0,
2189                                   Form("%s/HGLG_Corr_Layer%02d.%s" ,outputDirPlots.Data(), l, plotSuffix.Data()), it->second);
2190       }
2191     }
2192     std::cout << "done plotting single layers" << std::endl;  
2193   }
2194   return true;
2195 }
2196 
2197 //***********************************************************************************************
2198 //*********************** improved scaling calculation function *********************************
2199 //***********************************************************************************************
2200 bool Analyses::GetImprovedScaling(void){
2201   std::cout<<"GetImprovedScaling"<<std::endl;
2202   
2203   std::map<int,RunInfo> ri=readRunInfosFromFile(RunListInputName.Data(),debug,0);
2204   
2205   std::map<int,TileSpectra> hSpectra;
2206   std::map<int,TileSpectra> hSpectraTrigg;
2207   std::map<int, TileSpectra>::iterator ithSpectra;
2208   std::map<int, TileSpectra>::iterator ithSpectraTrigg;
2209   
2210   std::cout << "Additional Output with histos being created: " << RootOutputNameHist.Data() << std::endl;
2211   if(Overwrite){
2212     std::cout << "recreating file with hists" << std::endl;
2213     RootOutputHist = new TFile(RootOutputNameHist.Data(),"RECREATE");
2214   } else{
2215     std::cout << "newly creating file with hists" << std::endl;
2216     RootOutputHist = new TFile(RootOutputNameHist.Data(),"CREATE");
2217   }
2218     
2219   ROOT::Math::MinimizerOptions::SetDefaultMinimizer("Minuit2", "Migrad");
2220   // setup trigger sel
2221   double factorMinTrigg   = 0.8;
2222   double factorMaxTrigg   = 2.5;
2223   if (yearData == 2023){
2224     factorMinTrigg    = 0.9;
2225     factorMaxTrigg    = 2.;
2226   }
2227   
2228   RootOutputHist->mkdir("IndividualCells");
2229   RootOutputHist->mkdir("IndividualCellsTrigg");
2230   RootOutputHist->cd("IndividualCellsTrigg");  
2231   //***********************************************************************************************
2232   //************************* first pass over tree to extract spectra *****************************
2233   //***********************************************************************************************  
2234   TcalibIn->GetEntry(0);
2235   // check whether calib should be overwritten based on external text file
2236   if (OverWriteCalib){
2237     calib.ReadCalibFromTextFile(ExternalCalibFile,debug);
2238   }
2239   
2240   int evts=TdataIn->GetEntries();
2241   int runNr = -1;
2242   int actChI  = 0;
2243   double averageScale     = calib.GetAverageScaleHigh(actChI);
2244   double averageScaleLow  = calib.GetAverageScaleLow();
2245   std::cout << "average HG mip: " << averageScale << " LG mip: "<< averageScaleLow << "\t act. ch: "<< actChI << std::endl;
2246   for(int i=0; i<evts; i++){
2247     TdataIn->GetEntry(i);
2248     if (i == 0)runNr = event.GetRunNumber();
2249     TdataIn->GetEntry(i);
2250     if (i%5000 == 0 && i > 0 && debug > 0) std::cout << "Reading " <<  i << " / " << evts << " events" << std::endl;
2251     for(int j=0; j<event.GetNTiles(); j++){
2252       Caen* aTile=(Caen*)event.GetTile(j);
2253       if (i == 0 && debug > 2) std::cout << ((TString)setup->DecodeCellID(aTile->GetCellID())).Data() << std::endl;
2254       long currCellID = aTile->GetCellID();
2255       
2256       // read current tile
2257       ithSpectraTrigg=hSpectraTrigg.find(aTile->GetCellID());
2258       double hgCorr = aTile->GetADCHigh()-calib.GetPedestalMeanH(aTile->GetCellID());
2259       double lgCorr = aTile->GetADCLow()-calib.GetPedestalMeanL(aTile->GetCellID());
2260 
2261       // estimate local muon trigger
2262       bool localMuonTrigg = event.InspectIfLocalMuonTrigg(currCellID, averageScale, factorMinTrigg, factorMaxTrigg);
2263       
2264       if(ithSpectraTrigg!=hSpectraTrigg.end()){
2265         ithSpectraTrigg->second.FillTrigger(aTile->GetLocalTriggerPrimitive());
2266       } else {
2267         RootOutputHist->cd("IndividualCellsTrigg");
2268         hSpectraTrigg[currCellID]=TileSpectra("mipTrigg",currCellID,calib.GetTileCalib(currCellID),event.GetROtype(),debug);
2269         hSpectraTrigg[currCellID].FillTrigger(aTile->GetLocalTriggerPrimitive());;
2270         RootOutput->cd();
2271       }
2272       
2273       ithSpectra=hSpectra.find(aTile->GetCellID());
2274       if (ithSpectra!=hSpectra.end()){
2275         ithSpectra->second.FillSpectra(lgCorr,hgCorr);
2276         if (hgCorr > 3*calib.GetPedestalSigH(currCellID) && lgCorr > 3*calib.GetPedestalSigL(currCellID) && hgCorr < 3900 )
2277           ithSpectra->second.FillCorr(lgCorr,hgCorr);
2278       } else {
2279         RootOutputHist->cd("IndividualCells");
2280         hSpectra[currCellID]=TileSpectra("mip1st",currCellID,calib.GetTileCalib(currCellID),event.GetROtype(),debug);
2281         hSpectra[aTile->GetCellID()].FillSpectra(lgCorr,hgCorr);;
2282         if (hgCorr > 3*calib.GetPedestalSigH(aTile->GetCellID()) && lgCorr > 3*calib.GetPedestalSigL(aTile->GetCellID() && hgCorr < 3900) )
2283           hSpectra[aTile->GetCellID()].FillCorr(lgCorr,hgCorr);;
2284 
2285         RootOutput->cd();
2286       }
2287       
2288      
2289       // only fill tile spectra if 4 surrounding tiles on average are compatible with muon response
2290       if (localMuonTrigg){
2291         aTile->SetLocalTriggerBit(1);
2292         ithSpectraTrigg=hSpectraTrigg.find(aTile->GetCellID());
2293         ithSpectraTrigg->second.FillSpectra(lgCorr,hgCorr);
2294         if (hgCorr > 3*calib.GetPedestalSigH(currCellID) && lgCorr > 3*calib.GetPedestalSigL(currCellID) && hgCorr < 3900 )
2295           ithSpectraTrigg->second.FillCorr(lgCorr,hgCorr);
2296       }
2297     }
2298     RootOutput->cd();
2299     TdataOut->Fill();
2300   }
2301   TdataOut->Write();
2302   TsetupIn->CloneTree()->Write();
2303   
2304   //***********************************************************************************************
2305   //***** Monitoring histos for fits results of 2nd iteration ******************
2306   //***********************************************************************************************
2307   RootOutputHist->cd();
2308   int maxChannelPerLayer        = (setup->GetNMaxColumn()+1)*(setup->GetNMaxRow()+1);
2309   // monitoring trigger 
2310   TH2D* hmipTriggers              = new TH2D( "hmipTriggers","muon triggers; layer; brd channel; counts ",
2311                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
2312   hmipTriggers->SetDirectory(0);
2313   TH2D* hSuppresionNoise          = new TH2D( "hSuppresionNoise","S/B noise region; layer; brd channel; S/B noise region",
2314                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
2315   hSuppresionNoise->SetDirectory(0);
2316   TH2D* hSuppresionSignal         = new TH2D( "hSuppresionSignal","S/B signal region; layer; brd channel; S/B signal region",
2317                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
2318   hSuppresionSignal->SetDirectory(0);
2319 
2320   // monitoring 2nd iteration mip fits
2321   TH2D* hspectraHGMaxVsLayer   = new TH2D( "hspectraHGMaxVsLayer","Max High Gain; layer; brd channel; Max_{HG, 2^{nd}} (arb. units) ",
2322                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
2323   hspectraHGMaxVsLayer->SetDirectory(0);
2324   TH2D* hspectraHGFWHMVsLayer   = new TH2D( "hspectraHGFWHMVsLayer","FWHM High Gain; layer; brd channel; FWHM_{HG, 2^{nd}} (arb. units) ",
2325                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
2326   hspectraHGFWHMVsLayer->SetDirectory(0);
2327   TH2D* hspectraHGLMPVVsLayer   = new TH2D( "hspectraHGLMPVVsLayer","MPV High Gain; layer; brd channel; MPV_{HG, 2^{nd}} (arb. units) ",
2328                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
2329   hspectraHGLMPVVsLayer->SetDirectory(0);
2330   TH2D* hspectraHGLSigmaVsLayer = new TH2D( "hspectraHGLSigmaVsLayer","Sigma Landau High Gain; layer; brd channel; #sigma_{L,HG, 2^{nd}} (arb. units) ",
2331                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
2332   hspectraHGLSigmaVsLayer->SetDirectory(0);
2333   TH2D* hspectraHGGSigmaVsLayer = new TH2D( "hspectraHGGSigmaVsLayer","Sigma Gauss High Gain; layer; brd channel; #sigma_{G,HG, 2^{nd}} (arb. units) ",
2334                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
2335   hspectraHGGSigmaVsLayer->SetDirectory(0);
2336   TH2D* hspectraLGMaxVsLayer   = new TH2D( "hspectraLGMaxVsLayer","Max High Gain; layer; brd channel; Max_{LG, 2^{nd}} (arb. units) ",
2337                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
2338   hspectraLGMaxVsLayer->SetDirectory(0);
2339   TH2D* hspectraLGFWHMVsLayer   = new TH2D( "hspectraLGFWHMVsLayer","FWHM High Gain; layer; brd channel; FWHM_{LG, 2^{nd}} (arb. units) ",
2340                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
2341   hspectraLGFWHMVsLayer->SetDirectory(0);
2342   TH2D* hspectraLGLMPVVsLayer   = new TH2D( "hspectraLGLMPVVsLayer","MPV High Gain; layer; brd channel; MPV_{LG, 2^{nd}} (arb. units) ",
2343                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
2344   hspectraLGLMPVVsLayer->SetDirectory(0);
2345   TH2D* hspectraLGLSigmaVsLayer = new TH2D( "hspectraLGLSigmaVsLayer","Sigma Landau High Gain; layer; brd channel; #sigma_{L,LG, 2^{nd}} (arb. units) ",
2346                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
2347   hspectraLGLSigmaVsLayer->SetDirectory(0);
2348   TH2D* hspectraLGGSigmaVsLayer = new TH2D( "hspectraLGGSigmaVsLayer","Sigma Gauss High Gain; layer; brd channel; #sigma_{G,LG, 2^{nd}} (arb. units) ",
2349                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
2350   hspectraLGGSigmaVsLayer->SetDirectory(0);
2351 
2352   TH1D* hMaxHG             = new TH1D( "hMaxHG","Max High Gain ;Max_{HG} (arb. units) ; counts ",
2353                                             2000, -0.5, 2000-0.5);
2354   hMaxHG->SetDirectory(0);
2355   TH1D* hMaxLG             = new TH1D( "hMaxLG","Max Low Gain ;Max_{LG} (arb. units) ; counts ",
2356                                             400, -0.5, 400-0.5);
2357   hMaxLG->SetDirectory(0);
2358 
2359 
2360   int currCells = 0;
2361   double* parameters    = new double[6];
2362   double* parErrAndRes  = new double[6];
2363   bool isGood;
2364   double meanSB_NoiseR  = 0;
2365   double meanSB_SigR    = 0;
2366   if ( debug > 0)
2367     std::cout << "============================== start fitting improved iteration" << std::endl;  
2368   
2369   for(ithSpectraTrigg=hSpectraTrigg.begin(); ithSpectraTrigg!=hSpectraTrigg.end(); ++ithSpectraTrigg){
2370     if (currCells%20 == 0 && currCells > 0 && debug > 0)
2371       std::cout << "============================== cell " <<  currCells << " / " << hSpectraTrigg.size() << " cells" << std::endl;
2372     currCells++;
2373     for (int p = 0; p < 6; p++){
2374       parameters[p]   = 0;
2375       parErrAndRes[p] = 0;
2376     }
2377     isGood=ithSpectraTrigg->second.FitMipHG(parameters, parErrAndRes, debug, yearData, true, calib.GetVov(), averageScale);
2378     
2379     long cellID     = ithSpectraTrigg->second.GetCellID();
2380     int layer       = setup->GetLayer(cellID);
2381     int chInLayer   = setup->GetChannelInLayer(cellID);    
2382     int bin2D       = hspectraHGMaxVsLayer->FindBin(layer,chInLayer);
2383 
2384     Int_t binNLow   = ithSpectraTrigg->second.GetHG()->FindBin(-1*calib.GetPedestalSigH(cellID));
2385     Int_t binNHigh  = ithSpectraTrigg->second.GetHG()->FindBin(3*calib.GetPedestalSigH(cellID));
2386     Int_t binSHigh  = ithSpectraTrigg->second.GetHG()->FindBin(3800);
2387     
2388     double S_NoiseR = ithSpectraTrigg->second.GetHG()->Integral(binNLow, binNHigh);
2389     double S_SigR   = ithSpectraTrigg->second.GetHG()->Integral(binNHigh, binSHigh);
2390     
2391     ithSpectra      = hSpectra.find(cellID);
2392     double B_NoiseR = ithSpectra->second.GetHG()->Integral(binNLow , binNHigh);
2393     double B_SigR   = ithSpectra->second.GetHG()->Integral(binNHigh, binSHigh);
2394     
2395     double SB_NoiseR  = (B_NoiseR != 0.) ? S_NoiseR/B_NoiseR : 0;
2396     double SB_SigR    = (B_SigR != 0.) ? S_SigR/B_SigR : 0;
2397     
2398     meanSB_NoiseR += SB_NoiseR;
2399     meanSB_SigR += SB_SigR;
2400     
2401     hmipTriggers->SetBinContent(bin2D, ithSpectraTrigg->second.GetHG()->GetEntries());
2402     hSuppresionNoise->SetBinContent(bin2D, SB_NoiseR);
2403     hSuppresionSignal->SetBinContent(bin2D, SB_SigR);
2404     if (isGood){
2405       hspectraHGMaxVsLayer->SetBinContent(bin2D, parameters[4]);
2406       hspectraHGFWHMVsLayer->SetBinContent(bin2D, parameters[5]);
2407       hspectraHGLMPVVsLayer->SetBinContent(bin2D, parameters[1]);
2408       hspectraHGLMPVVsLayer->SetBinError(bin2D, parErrAndRes[1]);
2409       hspectraHGLSigmaVsLayer->SetBinContent(bin2D, parameters[0]);
2410       hspectraHGLSigmaVsLayer->SetBinError(bin2D, parErrAndRes[0]);
2411       hspectraHGGSigmaVsLayer->SetBinContent(bin2D, parameters[3]);
2412       hspectraHGGSigmaVsLayer->SetBinError(bin2D, parErrAndRes[3]);
2413       hMaxHG->Fill(parameters[4]);
2414     }
2415     for (int p = 0; p < 6; p++){
2416       parameters[p]   = 0;
2417       parErrAndRes[p] = 0;
2418     }
2419     isGood=ithSpectraTrigg->second.FitMipLG(parameters, parErrAndRes, debug, yearData, true, averageScaleLow);
2420     if (isGood){
2421       hspectraLGMaxVsLayer->SetBinContent(bin2D, parameters[4]);
2422       hspectraLGFWHMVsLayer->SetBinContent(bin2D, parameters[5]);
2423       hspectraLGLMPVVsLayer->SetBinContent(bin2D, parameters[1]);
2424       hspectraLGLMPVVsLayer->SetBinError(bin2D, parErrAndRes[1]);
2425       hspectraLGLSigmaVsLayer->SetBinContent(bin2D, parameters[0]);
2426       hspectraLGLSigmaVsLayer->SetBinError(bin2D, parErrAndRes[0]);
2427       hspectraLGGSigmaVsLayer->SetBinContent(bin2D, parameters[3]);
2428       hspectraLGGSigmaVsLayer->SetBinError(bin2D, parErrAndRes[3]);
2429       hMaxLG->Fill(parameters[4]);
2430     }
2431   }
2432   if ( debug > 0)
2433     std::cout << "============================== done fitting improved iteration" << std::endl;
2434 
2435   
2436   meanSB_NoiseR = meanSB_NoiseR/hSpectraTrigg.size();
2437   meanSB_SigR   = meanSB_SigR/hSpectraTrigg.size();
2438   
2439   RootOutput->cd();
2440   
2441   if (IsCalibSaveToFile()){
2442     TString fileCalibPrint = RootOutputName;
2443     fileCalibPrint         = fileCalibPrint.ReplaceAll(".root","_calib.txt");
2444     calib.PrintCalibToFile(fileCalibPrint);
2445   }
2446   
2447   TcalibOut->Fill();
2448   TcalibOut->Write();
2449   int actChA                     = 0;
2450   double averageScaleUpdated     = calib.GetAverageScaleHigh(actChA);
2451   double averageScaleUpdatedLow  = calib.GetAverageScaleLow();
2452   double meanFWHMHG              = calib.GetAverageScaleWidthHigh();
2453   double meanFWHMLG              = calib.GetAverageScaleWidthLow();
2454 
2455   std::cout << "average input HG mip: " << averageScale << " LG mip: "<< averageScaleLow << "\t act. ch: "<< actChI<< std::endl;
2456   std::cout << "average updated HG mip: " << averageScaleUpdated << " LG mip: "<< averageScaleUpdatedLow << "\t act. ch: "<< actChA<< std::endl;
2457 
2458   RootOutput->Close();
2459 
2460 
2461   RootOutputHist->cd("IndividualCellsTrigg");
2462     for(ithSpectra=hSpectraTrigg.begin(); ithSpectra!=hSpectraTrigg.end(); ++ithSpectra){
2463       ithSpectra->second.Write(true);
2464     }
2465   RootOutputHist->cd();
2466     
2467     hspectraHGMaxVsLayer->Write();
2468     hspectraHGFWHMVsLayer->Write();
2469     hspectraHGLMPVVsLayer->Write();
2470     hspectraHGLSigmaVsLayer->Write();
2471     hspectraHGGSigmaVsLayer->Write();
2472     hMaxHG->Write();
2473     
2474     hspectraLGMaxVsLayer->Write();
2475     hspectraLGFWHMVsLayer->Write();
2476     hspectraLGLMPVVsLayer->Write();
2477     hspectraLGLSigmaVsLayer->Write();
2478     hspectraLGGSigmaVsLayer->Write();
2479     hMaxLG->Write();
2480     
2481     hmipTriggers->Write();
2482     hSuppresionNoise->Write();
2483     hSuppresionSignal->Write();
2484   // fill calib tree & write it
2485   // close open root files
2486   RootOutputHist->Write();
2487   RootOutputHist->Close();
2488 
2489   RootInput->Close();
2490 
2491   // Get run info object
2492   std::map<int,RunInfo>::iterator it=ri.find(runNr);
2493   
2494   // create directory for plot output
2495   TString outputDirPlots = GetPlotOutputDir();
2496   gSystem->Exec("mkdir -p "+outputDirPlots);
2497   
2498   //**********************************************************************
2499   // Create canvases for channel overview plotting
2500   //**********************************************************************
2501   Double_t textSizeRel = 0.035;
2502   StyleSettingsBasics("pdf");
2503   SetPlotStyle();
2504 
2505   TCanvas* canvas2DCorr = new TCanvas("canvasCorrPlots","",0,0,1450,1300);  // gives the page size
2506   DefaultCancasSettings( canvas2DCorr, 0.08, 0.13, 0.045, 0.07);
2507 
2508   canvas2DCorr->SetLogz(0);
2509   PlotSimple2D( canvas2DCorr, hspectraHGMaxVsLayer, -10000, -10000, textSizeRel, Form("%s/HG_MaxMip.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true, Form( "#LT Max_{HG} #GT = %.1f", averageScaleUpdated) );
2510   PlotSimple2D( canvas2DCorr, hspectraHGFWHMVsLayer, -10000, -10000, textSizeRel, Form("%s/HG_FWHMMip.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true, Form( "#LT FWHM_{HG} #GT = %.1f", meanFWHMHG));
2511   PlotSimple2D( canvas2DCorr, hspectraHGLMPVVsLayer, -10000, -10000, textSizeRel, Form("%s/HG_LandMPVMip.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
2512   PlotSimple2D( canvas2DCorr, hspectraHGLSigmaVsLayer, -10000, -10000, textSizeRel, Form("%s/HG_LandSigMip.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
2513   PlotSimple2D( canvas2DCorr, hspectraHGGSigmaVsLayer, -10000, -10000, textSizeRel, Form("%s/HG_GaussSigMip.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
2514   canvas2DCorr->SetLogz(1);
2515   TString drawOpt = "colz";
2516   if (yearData == 2023)
2517     drawOpt = "colz,text";
2518   PlotSimple2D( canvas2DCorr, hmipTriggers, -10000, -10000, textSizeRel, Form("%s/MuonTriggers.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, drawOpt, true, Form( "evt. collected = %d", evts));
2519   PlotSimple2D( canvas2DCorr, hSuppresionNoise, -10000, -10000, textSizeRel, Form("%s/SuppressionNoise.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, drawOpt, true, Form( "#LT S/B noise #GT = %.3f", meanSB_NoiseR));
2520   PlotSimple2D( canvas2DCorr, hSuppresionSignal, -10000, -10000, textSizeRel, Form("%s/SuppressionSignal.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, drawOpt, true, Form( "#LT S/B signal #GT = %.3f", meanSB_SigR));
2521   
2522   canvas2DCorr->SetLogz(0);
2523   PlotSimple2D( canvas2DCorr, hspectraLGMaxVsLayer, -10000, -10000, textSizeRel, Form("%s/LG_MaxMip.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true, Form( "#LT Max_{LG} #GT = %.1f", averageScaleUpdatedLow));
2524   PlotSimple2D( canvas2DCorr, hspectraLGFWHMVsLayer, -10000, -10000, textSizeRel, Form("%s/LG_FWHMMip.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true, Form( "#LT FWHM_{LG} #GT = %.1f", meanFWHMLG));
2525   PlotSimple2D( canvas2DCorr, hspectraLGLMPVVsLayer, -10000, -10000, textSizeRel, Form("%s/LG_LandMPVMip.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
2526   PlotSimple2D( canvas2DCorr, hspectraLGLSigmaVsLayer, -10000, -10000, textSizeRel, Form("%s/LG_LandSigMip.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
2527   PlotSimple2D( canvas2DCorr, hspectraLGGSigmaVsLayer, -10000, -10000, textSizeRel, Form("%s/LG_GaussSigMip.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
2528 
2529   //***********************************************************************************************************
2530   //********************************* 8 Panel overview plot  **************************************************
2531   //***********************************************************************************************************
2532   //*****************************************************************
2533     // Test beam geometry (beam coming from viewer)
2534     //===========================================================
2535     //||    8 (4)    ||    7 (5)   ||    6 (6)   ||    5 (7)   ||  row 0
2536     //===========================================================
2537     //||    1 (0)    ||    2 (1)   ||    3 (2)   ||    4 (3)   ||  row 1
2538     //===========================================================
2539     //    col 0     col 1       col 2     col  3
2540     // rebuild pad geom in similar way (numbering -1)
2541   //*****************************************************************
2542   TCanvas* canvas8Panel;
2543   TPad* pad8Panel[8];
2544   Double_t topRCornerX[8];
2545   Double_t topRCornerY[8];
2546   Int_t textSizePixel = 30;
2547   Double_t relSize8P[8];
2548   CreateCanvasAndPadsFor8PannelTBPlot(canvas8Panel, pad8Panel,  topRCornerX, topRCornerY, relSize8P, textSizePixel);
2549  
2550   calib.PrintGlobalInfo();  
2551   Double_t maxHG = ReturnMipPlotRangeDepVov(calib.GetVov(),true);
2552   Double_t maxLG = ReturnMipPlotRangeDepVov(calib.GetVov(),false);
2553   std::cout << "plotting single layers" << std::endl;
2554   for (Int_t l = 0; l < setup->GetNMaxLayer()+1; l++){    
2555     if (l%10 == 0 && l > 0 && debug > 0)
2556       std::cout << "============================== layer " <<  l << " / " << setup->GetNMaxLayer() << " layers" << std::endl;
2557     PlotMipWithFitsFullLayer (canvas8Panel,pad8Panel, topRCornerX, topRCornerY, relSize8P, textSizePixel, 
2558                               hSpectra, hSpectraTrigg, setup, true, -100, maxHG, 1.2, l, 0,
2559                               Form("%s/MIP_HG_Layer%02d.%s" ,outputDirPlots.Data(), l, plotSuffix.Data()), it->second);
2560     PlotMipWithFitsFullLayer (canvas8Panel,pad8Panel, topRCornerX, topRCornerY, relSize8P, textSizePixel, 
2561                               hSpectra, hSpectraTrigg, setup, false, -20, maxLG, 1.2, l, 0,
2562                               Form("%s/MIP_LG_Layer%02d.%s" ,outputDirPlots.Data(), l, plotSuffix.Data()), it->second);
2563     PlotTriggerPrimWithFitsFullLayer (canvas8Panel,pad8Panel, topRCornerX, topRCornerY, relSize8P, textSizePixel, 
2564                                       hSpectraTrigg, setup, averageScale, factorMinTrigg, factorMaxTrigg,
2565                                       0, maxHG*2, 1.2, l, 0, Form("%s/TriggPrimitive_Layer%02d.%s" ,outputDirPlots.Data(), l, plotSuffix.Data()), it->second);
2566   }
2567   std::cout << "done plotting" << std::endl;
2568   
2569   
2570   return true;
2571 }
2572 
2573 
2574 
2575 //***********************************************************************************************
2576 //*********************** improved scaling calculation function *********************************
2577 //***********************************************************************************************
2578 bool Analyses::GetNoiseSampleAndRefitPedestal(void){
2579   std::cout<<"GetNoiseSampleAndRefitPedestal"<<std::endl;
2580   
2581   std::map<int,RunInfo> ri=readRunInfosFromFile(RunListInputName.Data(),debug,0);
2582   
2583   std::map<int,TileSpectra> hSpectra;
2584   std::map<int,TileSpectra> hSpectraTrigg;
2585   std::map<int, TileSpectra>::iterator ithSpectra;
2586   std::map<int, TileSpectra>::iterator ithSpectraTrigg;
2587   
2588   std::cout << "Additional Output with histos being created: " << RootOutputNameHist.Data() << std::endl;
2589   if(Overwrite){
2590     std::cout << "recreating file with hists" << std::endl;
2591     RootOutputHist = new TFile(RootOutputNameHist.Data(),"RECREATE");
2592   } else{
2593     std::cout << "newly creating file with hists" << std::endl;
2594     RootOutputHist = new TFile(RootOutputNameHist.Data(),"CREATE");
2595   }
2596     
2597   // setup trigger sel
2598   double factorMinTrigg   = 0.5;
2599   if(yearData == 2023)
2600     factorMinTrigg        = 0.1;
2601   // create HG and LG histo's per channel
2602   TH2D* hspectraHGvsCellID      = new TH2D( "hNoiseTriggeredSpectraHG_vsCellID","Noise trigg ADC spectrum High Gain vs CellID; cell ID; ADC_{HG} (arb. units) ",
2603                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5, 4000,0,4000);
2604   hspectraHGvsCellID->SetDirectory(0);
2605   TH2D* hspectraLGvsCellID      = new TH2D( "hNoiseTriggeredSpectraLG_vsCellID","Noise trigg ADC spectrum Low Gain vs CellID; cell ID; ADC_{LG} (arb. units)  ",
2606                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5, 4000,0,4000);
2607   hspectraLGvsCellID->SetDirectory(0);
2608 
2609   
2610   RootOutputHist->mkdir("IndividualCells");
2611   RootOutputHist->mkdir("IndividualCellsTrigg");
2612   RootOutputHist->cd("IndividualCellsTrigg");  
2613   
2614   //***********************************************************************************************
2615   //************************* first pass over tree to extract spectra *****************************
2616   //***********************************************************************************************  
2617   TcalibIn->GetEntry(0);
2618   // check whether calib should be overwritten based on external text file
2619   if (OverWriteCalib){
2620     calib.ReadCalibFromTextFile(ExternalCalibFile,debug);
2621   }
2622   
2623   int evts=TdataIn->GetEntries();
2624   int runNr = -1;
2625   int actCh = 0;
2626   double averageScale     = calib.GetAverageScaleHigh(actCh);
2627   double averageScaleLow  = calib.GetAverageScaleLow();
2628   std::cout << "average HG mip: " << averageScale << " LG mip: "<< averageScaleLow << std::endl;
2629   for(int i=0; i<evts; i++){
2630     TdataIn->GetEntry(i);
2631     if (i == 0)runNr = event.GetRunNumber();
2632     TdataIn->GetEntry(i);
2633     if (i%5000 == 0 && i > 0 && debug > 0) std::cout << "Reading " <<  i << " / " << evts << " events" << std::endl;
2634     for(int j=0; j<event.GetNTiles(); j++){
2635       Caen* aTile=(Caen*)event.GetTile(j);
2636       if (i == 0 && debug > 2) std::cout << ((TString)setup->DecodeCellID(aTile->GetCellID())).Data() << std::endl;
2637       long currCellID = aTile->GetCellID();
2638       
2639       // read current tile
2640       ithSpectraTrigg=hSpectraTrigg.find(aTile->GetCellID());
2641       // estimate local muon trigger
2642       bool localNoiseTrigg = event.InspectIfNoiseTrigg(currCellID, averageScale, factorMinTrigg);
2643       
2644       if(ithSpectraTrigg!=hSpectraTrigg.end()){
2645         ithSpectraTrigg->second.FillTrigger(aTile->GetLocalTriggerPrimitive());
2646       } else {
2647         RootOutputHist->cd("IndividualCellsTrigg");
2648         hSpectraTrigg[currCellID]=TileSpectra("mipTrigg",currCellID,calib.GetTileCalib(currCellID),event.GetROtype(),debug);
2649         hSpectraTrigg[currCellID].FillTrigger(aTile->GetLocalTriggerPrimitive());;
2650         RootOutput->cd();
2651       }
2652       
2653       ithSpectra=hSpectra.find(aTile->GetCellID());
2654       if (ithSpectra!=hSpectra.end()){
2655         ithSpectra->second.FillSpectra(aTile->GetADCLow(),aTile->GetADCHigh());
2656       } else {
2657         RootOutputHist->cd("IndividualCells");
2658         hSpectra[currCellID]=TileSpectra("mip1st",currCellID,calib.GetTileCalib(currCellID),event.GetROtype(),debug);
2659         hSpectra[aTile->GetCellID()].FillSpectra(aTile->GetADCLow(),aTile->GetADCHigh());;
2660 
2661         RootOutput->cd();
2662       }
2663      
2664       // only fill tile spectra if X surrounding tiles on average are compatible with pure noise
2665       if (localNoiseTrigg){
2666         aTile->SetLocalTriggerBit(2);
2667         ithSpectraTrigg=hSpectraTrigg.find(aTile->GetCellID());
2668         ithSpectraTrigg->second.FillSpectra(aTile->GetADCLow(),aTile->GetADCHigh());
2669         
2670         hspectraHGvsCellID->Fill(aTile->GetCellID(), aTile->GetADCHigh());
2671         hspectraLGvsCellID->Fill(aTile->GetCellID(), aTile->GetADCLow());
2672       }
2673     }
2674     RootOutput->cd();
2675     TdataOut->Fill();
2676   }
2677   TdataOut->Write();
2678   TsetupIn->CloneTree()->Write();
2679   
2680   //***********************************************************************************************
2681   //***** Monitoring histos for fits results of 2nd iteration ******************
2682   //***********************************************************************************************
2683   RootOutputHist->cd();
2684   int maxChannelPerLayer        = (setup->GetNMaxColumn()+1)*(setup->GetNMaxRow()+1);
2685   // monitoring trigger 
2686   TH2D* hnoiseTriggers              = new TH2D( "hnoiseTriggers","muon triggers; layer; brd channel; counts ",
2687                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
2688   hnoiseTriggers->SetDirectory(0);
2689   TH1D* hMeanPedHGvsCellID      = new TH1D( "hMeanPedHG_vsCellID","mean Ped High Gain vs CellID ; cell ID; #mu_{noise, HG} (arb. units) ",
2690                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5);
2691   hMeanPedHGvsCellID->SetDirectory(0);
2692   TH2D* hspectraHGMeanVsLayer   = new TH2D( "hspectraHGMeanVsLayer","Mean Ped High Gain vs CellID; layer; brd channel; #mu_{Ped HG} (arb. units) ",
2693                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
2694   hspectraHGMeanVsLayer->SetDirectory(0);
2695   TH2D* hspectraHGSigmaVsLayer  = new TH2D( "hspectraHGSigmaVsLayer","Mean Ped High Gain vs CellID; layer; brd channel; #sigma_{Ped HG} (arb. units) ",
2696                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
2697   hspectraHGSigmaVsLayer->SetDirectory(0);
2698   TH1D* hMeanPedLGvsCellID      = new TH1D( "hMeanPedLG_vsCellID","mean Ped Low Gain vs CellID ; cell ID; #mu_{noise, LG} (arb. units) ",
2699                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5);
2700   hMeanPedLGvsCellID->SetDirectory(0);
2701   TH2D* hspectraLGMeanVsLayer   = new TH2D( "hspectraLGMeanVsLayer","Mean Ped Low Gain vs CellID; layer; brd channel; #mu_{PED LG} (arb. units) ",
2702                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
2703   hspectraLGMeanVsLayer->SetDirectory(0);
2704   TH2D* hspectraLGSigmaVsLayer  = new TH2D( "hspectraLGSigmaVsLayer","Mean Ped Low Gain vs CellID; layer; brd channel; #sigma_{Ped LG} (arb. units)",
2705                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
2706   hspectraLGSigmaVsLayer->SetDirectory(0);
2707 
2708   if ( debug > 0)
2709     std::cout << "============================== starting fitting" << std::endl;
2710 
2711   int currCells = 0;
2712   double* parameters    = new double[6];
2713   for(ithSpectraTrigg=hSpectraTrigg.begin(); ithSpectraTrigg!=hSpectraTrigg.end(); ++ithSpectraTrigg){
2714     for (int p = 0; p < 6; p++){
2715       parameters[p]   = 0;
2716     }
2717 
2718     if (currCells%20 == 0 && currCells > 0 && debug > 0)
2719       std::cout << "============================== cell " <<  currCells << " / " << hSpectraTrigg.size() << " cells" << std::endl;
2720     currCells++;
2721     if ( debug > 2) std::cout << ((TString)setup->DecodeCellID(ithSpectraTrigg->second.GetCellID())).Data() << std::endl;
2722     ithSpectraTrigg->second.FitNoise(parameters, yearData, true);
2723     hMeanPedHGvsCellID->SetBinContent(hMeanPedHGvsCellID->FindBin(ithSpectraTrigg->second.GetCellID()), parameters[4]);
2724     hMeanPedHGvsCellID->SetBinError  (hMeanPedHGvsCellID->FindBin(ithSpectraTrigg->second.GetCellID()), parameters[6]);
2725     hMeanPedLGvsCellID->SetBinContent(hMeanPedLGvsCellID->FindBin(ithSpectraTrigg->second.GetCellID()), parameters[0]);
2726     hMeanPedLGvsCellID->SetBinError  (hMeanPedLGvsCellID->FindBin(ithSpectraTrigg->second.GetCellID()), parameters[2]);
2727     
2728     int layer     = setup->GetLayer(ithSpectraTrigg->second.GetCellID());
2729     int chInLayer = setup->GetChannelInLayer(ithSpectraTrigg->second.GetCellID());
2730   
2731     hspectraHGMeanVsLayer->SetBinContent(hspectraHGMeanVsLayer->FindBin(layer,chInLayer), parameters[4]);
2732     hspectraHGMeanVsLayer->SetBinError(hspectraHGMeanVsLayer->FindBin(layer,chInLayer), parameters[5]);
2733     hspectraHGSigmaVsLayer->SetBinContent(hspectraHGSigmaVsLayer->FindBin(layer,chInLayer), parameters[6]);
2734     hspectraHGSigmaVsLayer->SetBinError(hspectraHGSigmaVsLayer->FindBin(layer,chInLayer), parameters[7]);
2735     hspectraLGMeanVsLayer->SetBinContent(hspectraLGMeanVsLayer->FindBin(layer,chInLayer), parameters[0]);
2736     hspectraLGMeanVsLayer->SetBinError(hspectraLGMeanVsLayer->FindBin(layer,chInLayer), parameters[1]);
2737     hspectraLGSigmaVsLayer->SetBinContent(hspectraLGSigmaVsLayer->FindBin(layer,chInLayer), parameters[2]);
2738     hspectraLGSigmaVsLayer->SetBinError(hspectraLGSigmaVsLayer->FindBin(layer,chInLayer), parameters[3]);
2739     
2740     hnoiseTriggers->SetBinContent(hnoiseTriggers->FindBin(layer,chInLayer), ithSpectraTrigg->second.GetHG()->GetEntries());
2741   }
2742   if ( debug > 0)
2743     std::cout << "============================== done fitting" << std::endl;
2744   
2745   RootOutput->cd();
2746   
2747   if (IsCalibSaveToFile()){
2748     TString fileCalibPrint = RootOutputName;
2749     fileCalibPrint         = fileCalibPrint.ReplaceAll(".root","_calib.txt");
2750     calib.PrintCalibToFile(fileCalibPrint);
2751   }
2752 
2753 
2754   TcalibOut->Fill();
2755   TcalibOut->Write();
2756   
2757   RootOutput->Write();
2758   RootOutput->Close();
2759 
2760 
2761   RootOutputHist->cd("IndividualCellsTrigg");
2762     for(ithSpectra=hSpectraTrigg.begin(); ithSpectra!=hSpectraTrigg.end(); ++ithSpectra){
2763       ithSpectra->second.Write(true);
2764     }
2765   RootOutputHist->cd();
2766     
2767     hMeanPedHGvsCellID->Write();
2768     hMeanPedLGvsCellID->Write();
2769     hspectraHGMeanVsLayer->Write();
2770     hspectraHGSigmaVsLayer->Write();
2771     hspectraLGMeanVsLayer->Write(); 
2772     hspectraLGSigmaVsLayer->Write();
2773     hspectraHGvsCellID->Write();
2774     hspectraLGvsCellID->Write();
2775         
2776     hnoiseTriggers->Write();
2777   // close open root files
2778   RootOutputHist->Write();
2779   RootOutputHist->Close();
2780 
2781   RootInput->Close();
2782 
2783   // Get run info object
2784   std::map<int,RunInfo>::iterator it=ri.find(runNr);
2785   
2786   // create directory for plot output
2787   TString outputDirPlots = GetPlotOutputDir();
2788   gSystem->Exec("mkdir -p "+outputDirPlots);
2789   
2790   //**********************************************************************
2791   // Create canvases for channel overview plotting
2792   //**********************************************************************
2793   Double_t textSizeRel = 0.035;
2794   StyleSettingsBasics("pdf");
2795   SetPlotStyle();
2796 
2797   TCanvas* canvas2DCorr = new TCanvas("canvasCorrPlots","",0,0,1450,1300);  // gives the page size
2798   DefaultCancasSettings( canvas2DCorr, 0.08, 0.13, 0.045, 0.07);
2799 
2800   canvas2DCorr->SetLogz(0);
2801   PlotSimple2D( canvas2DCorr, hspectraHGMeanVsLayer, -10000, -10000, textSizeRel, Form("%s/HG_NoiseMean.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true );
2802   PlotSimple2D( canvas2DCorr, hspectraHGSigmaVsLayer, -10000, -10000, textSizeRel, Form("%s/HG_NoiseSigma.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
2803   PlotSimple2D( canvas2DCorr, hspectraLGMeanVsLayer, -10000, -10000, textSizeRel, Form("%s/LG_NoiseMean.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
2804   PlotSimple2D( canvas2DCorr, hspectraLGSigmaVsLayer, -10000, -10000, textSizeRel, Form("%s/LG_NoiseSigma.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
2805   canvas2DCorr->SetLogz(1);
2806   PlotSimple2D( canvas2DCorr, hspectraHGvsCellID, -10000, -10000, textSizeRel, Form("%s/HG_Noise.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
2807   PlotSimple2D( canvas2DCorr, hspectraLGvsCellID, -10000, -10000, textSizeRel, Form("%s/LG_Noise.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
2808   
2809   PlotSimple2D( canvas2DCorr, hnoiseTriggers, -10000, -10000, textSizeRel, Form("%s/LG_Noise.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true, Form( "evt. coll = %d", evts));
2810   //***********************************************************************************************************
2811   //********************************* 8 Panel overview plot  **************************************************
2812   //***********************************************************************************************************
2813   //*****************************************************************
2814     // Test beam geometry (beam coming from viewer)
2815     //===========================================================
2816     //||    8 (4)    ||    7 (5)   ||    6 (6)   ||    5 (7)   ||  row 0
2817     //===========================================================
2818     //||    1 (0)    ||    2 (1)   ||    3 (2)   ||    4 (3)   ||  row 1
2819     //===========================================================
2820     //    col 0     col 1       col 2     col  3
2821     // rebuild pad geom in similar way (numbering -1)
2822   //*****************************************************************
2823   TCanvas* canvas8Panel;
2824   TPad* pad8Panel[8];
2825   Double_t topRCornerX[8];
2826   Double_t topRCornerY[8];
2827   Int_t textSizePixel = 30;
2828   Double_t relSize8P[8];
2829   CreateCanvasAndPadsFor8PannelTBPlot(canvas8Panel, pad8Panel,  topRCornerX, topRCornerY, relSize8P, textSizePixel);
2830  
2831   for (Int_t l = 0; l < setup->GetNMaxLayer()+1; l++){    
2832     PlotNoiseAdvWithFitsFullLayer (canvas8Panel,pad8Panel, topRCornerX, topRCornerY, relSize8P, textSizePixel, 
2833                                     hSpectra, hSpectraTrigg, true, 0, 450, 1.2, l, 0,
2834                                     Form("%s/NoiseTrigg_HG_Layer%02d.%s" ,outputDirPlots.Data(), l, plotSuffix.Data()), it->second);
2835   }
2836 
2837   
2838   return true;
2839 }
2840 
2841 
2842 //***********************************************************************************************
2843 //*********************** Calibration routine ***************************************************
2844 //***********************************************************************************************
2845 bool Analyses::RunEvalLocalTriggers(void){
2846   std::cout<<"EvalLocalTriggers"<<std::endl;
2847 
2848   RootOutput->cd();
2849   std::cout << "starting to run trigger eval: " << TcalibIn <<  "\t" << TcalibIn->GetEntry(0) << std::endl;
2850   TcalibIn->GetEntry(0);
2851   // check whether calib should be overwritten based on external text file
2852   if (OverWriteCalib){
2853     calib.ReadCalibFromTextFile(ExternalCalibFile,debug);
2854   }
2855   
2856   int actCh1st = 0;
2857   double averageScale = calib.GetAverageScaleHigh(actCh1st);
2858   double avLGHGCorr   = calib.GetAverageLGHGCorr();
2859   std::cout << "average HG mip: " << averageScale << "\t active ch: "<< actCh1st<< std::endl;
2860   
2861   // setup local trigger sel
2862   TRandom3* rand    = new TRandom3();
2863   Int_t localTriggerTiles = 4;
2864   if (yearData == 2023){
2865     localTriggerTiles = 6;
2866   }
2867   double factorMinTrigg   = 0.8;
2868   double factorMinTriggNoise = 0.2;
2869   double factorMaxTrigg   = 2.;
2870   if (yearData == 2023){
2871     factorMinTrigg    = 0.9;
2872     factorMaxTrigg    = 2.;
2873   }
2874   
2875   int outCount      = 1000;
2876   int evts=TdataIn->GetEntries();
2877   
2878   if (maxEvents == -1){
2879     maxEvents = evts;
2880   } else {
2881     std::cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl;
2882     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;
2883     std::cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl;
2884   }
2885   
2886   if (evts < 10000)
2887     outCount  = 500;
2888   if (evts > 100000)
2889     outCount  = 5000;
2890   int runNr = -1;
2891   for(int i=0; i<evts && i < maxEvents; i++){
2892     TdataIn->GetEntry(i);
2893     if (i == 0){
2894       runNr = event.GetRunNumber();
2895       std::cout<< "original run numbers calib: "<<calib.GetRunNumber() << "\t" << calib.GetRunNumberPed() << "\t" << calib.GetRunNumberMip() << std::endl;
2896       calib.SetRunNumber(runNr);
2897       calib.SetBeginRunTime(event.GetBeginRunTimeAlt());
2898       std::cout<< "reset run numbers calib: "<< calib.GetRunNumber() << "\t" << calib.GetRunNumberPed() << "\t" << calib.GetRunNumberMip() << std::endl;
2899     }
2900     
2901     if (i%outCount == 0 && debug > 0) std::cout << "Reading " <<  i << " / " << evts << " events" << std::endl;
2902     for(int j=0; j<event.GetNTiles(); j++){
2903       Caen* aTile=(Caen*)event.GetTile(j);      
2904       // calculate trigger primitives
2905       aTile->SetLocalTriggerPrimitive(event.CalculateLocalMuonTrigg(calib, rand, aTile->GetCellID(), localTriggerTiles, avLGHGCorr));
2906       bool localMuonTrigg   = event.InspectIfLocalMuonTrigg(aTile->GetCellID(), averageScale, factorMinTrigg, factorMaxTrigg);
2907       bool localNoiseTrigg  = event.InspectIfNoiseTrigg(aTile->GetCellID(), averageScale, factorMinTriggNoise);
2908       aTile->SetLocalTriggerBit(0);
2909       if (localMuonTrigg) aTile->SetLocalTriggerBit(1);
2910       if (localNoiseTrigg) aTile->SetLocalTriggerBit(2);
2911     }
2912     TdataOut->Fill();
2913   }
2914   TdataOut->Write();
2915   TsetupIn->CloneTree()->Write();
2916   
2917   if (IsCalibSaveToFile()){
2918     TString fileCalibPrint = RootOutputName;
2919     fileCalibPrint         = fileCalibPrint.ReplaceAll(".root","_calib.txt");
2920     calib.PrintCalibToFile(fileCalibPrint);
2921   }
2922   
2923   TcalibOut->Fill();
2924   TcalibOut->Write();
2925   
2926   RootOutput->Close();
2927   RootInput->Close();      
2928   
2929   return true;
2930 }
2931 
2932 //***********************************************************************************************
2933 //*********************** Calibration routine ***************************************************
2934 //***********************************************************************************************
2935 bool Analyses::Calibrate(void){
2936   std::cout<<"Calibrate"<<std::endl;
2937 
2938   std::map<int,RunInfo> ri=readRunInfosFromFile(RunListInputName.Data(),debug,0);
2939 
2940   std::cout << "Additional Output with histos being created: " << RootOutputNameHist.Data() << std::endl;
2941   if(Overwrite){
2942     std::cout << "recreating file with hists" << std::endl;
2943     RootOutputHist = new TFile(RootOutputNameHist.Data(),"RECREATE");
2944   } else{
2945     std::cout << "newly creating file with hists" << std::endl;
2946     RootOutputHist = new TFile(RootOutputNameHist.Data(),"CREATE");
2947   }
2948   
2949   // create HG and LG histo's per channel
2950   TH2D* hspectraHGCorrvsCellID      = new TH2D( "hspectraHGCorr_vsCellID","ADC spectrum High Gain corrected vs CellID; cell ID; ADC_{HG} (arb. units)  ; counts ",
2951                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5, 4000,-200,3800);
2952   hspectraHGCorrvsCellID->SetDirectory(0);
2953   TH2D* hspectraHGCorrvsCellIDNoise      = new TH2D( "hspectraHGCorr_vsCellID_Noise","ADC spectrum High Gain corrected vs CellID Noise; cell ID; ADC_{HG} (arb. units)  ; counts ",
2954                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5, 4000,-200,3800);
2955   hspectraHGCorrvsCellIDNoise->SetDirectory(0);
2956   TH2D* hspectraLGCorrvsCellID      = new TH2D( "hspectraLGCorr_vsCellID","ADC spectrum Low Gain corrected vs CellID; cell ID; ADC_{LG} (arb. units)  ; counts  ",
2957                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5, 4000,-200,3800);
2958   hspectraLGCorrvsCellID->SetDirectory(0);
2959   TH2D* hspectraLGCorrvsCellIDNoise      = new TH2D( "hspectraLGCorr_vsCellID_Noise","ADC spectrum Low Gain corrected vs CellID Noise; cell ID; ADC_{LG} (arb. units)  ; counts  ",
2960                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5, 4000,-200,3800);
2961   hspectraLGCorrvsCellIDNoise->SetDirectory(0);
2962   TH2D* hspectraHGvsCellID      = new TH2D( "hspectraHG_vsCellID","ADC spectrum High Gain vs CellID; cell ID; ADC_{HG} (arb. units)   ; counts ",
2963                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5, 4000,0,4000);
2964   hspectraHGvsCellID->SetDirectory(0);
2965   TH2D* hspectraLGvsCellID      = new TH2D( "hspectraLG_vsCellID","ADC spectrum Low Gain vs CellID; cell ID; ADC_{LG} (arb. units)  ; counts",
2966                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5, 4000,0,4000);
2967   hspectraLGvsCellID->SetDirectory(0);
2968   TH2D* hspectraEnergyvsCellID  = new TH2D( "hspectraEnergy_vsCellID","Energy vs CellID; cell ID; E (mip eq./tile)    ; counts",
2969                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5, 6000,0,200);
2970   hspectraEnergyvsCellID->SetDirectory(0);
2971   TH2D* hspectraEnergyTotvsNCells  = new TH2D( "hspectraTotEnergy_vsNCells","Energy vs CellID; N_{Cells}; E_{tot} (mip eq./tile) ; counts",
2972                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5, 6000,0,1000);
2973   hspectraEnergyTotvsNCells->SetDirectory(0);
2974 
2975   std::map<int,TileSpectra> hSpectra;
2976   std::map<int, TileSpectra>::iterator ithSpectra;
2977   std::map<int,TileSpectra> hSpectraNoise;
2978   std::map<int, TileSpectra>::iterator ithSpectraNoise;
2979   
2980   // entering histoOutput file
2981   RootOutputHist->mkdir("IndividualCells");
2982   RootOutputHist->cd("IndividualCells");
2983   RootOutputHist->mkdir("IndividualCellsNoise");
2984   RootOutputHist->cd("IndividualCellsNoise");
2985 
2986   Int_t runNr = -1;
2987   RootOutput->cd();
2988   std::cout << "starting to run calibration: " << TcalibIn <<  "\t" << TcalibIn->GetEntry(0) << std::endl;
2989   TcalibIn->GetEntry(0);
2990   // check whether calib should be overwritten based on external text file
2991   if (OverWriteCalib){
2992     calib.ReadCalibFromTextFile(ExternalCalibFile,debug);
2993   }
2994 
2995   int actCh1st = 0;
2996   double averageScale = calib.GetAverageScaleHigh(actCh1st);
2997   double avLGHGCorr   = calib.GetAverageLGHGCorr();
2998   std::cout << "average HG mip: " << averageScale << "\t active ch: "<< actCh1st<< std::endl;
2999   
3000   // setup local trigger sel
3001   TRandom3* rand    = new TRandom3();
3002   Int_t localTriggerTiles = 4;
3003   if (yearData == 2023){
3004     localTriggerTiles = 6;
3005   }
3006   double factorMinTrigg   = 0.8;
3007   double factorMinTriggNoise = 0.2;
3008   double factorMaxTrigg   = 2.;
3009   if (yearData == 2023){
3010     factorMinTrigg    = 0.9;
3011     factorMaxTrigg    = 2.;
3012   }
3013 
3014   
3015   double minMipFrac = 0.3;
3016   int corrHGADCSwap = 3500;
3017   int outCount      = 5000;
3018   int evts=TdataIn->GetEntries();
3019   
3020   if (maxEvents == -1){
3021     maxEvents = evts;
3022   } else {
3023     std::cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl;
3024     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;
3025     std::cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl;
3026   }
3027 
3028   if (evts < 10000)
3029     outCount  = 500;
3030   for(int i=0; i<evts && i < maxEvents; i++){
3031     TdataIn->GetEntry(i);
3032     if (i%outCount == 0 && debug > 0) std::cout << "Reading " <<  i << " / " << evts << " events" << std::endl;
3033     if (i == 0){
3034       runNr = event.GetRunNumber();
3035       std::cout<< "original run numbers calib: "<<calib.GetRunNumber() << "\t" << calib.GetRunNumberPed() << "\t" << calib.GetRunNumberMip() << std::endl;
3036       calib.SetRunNumber(runNr);
3037       calib.SetBeginRunTime(event.GetBeginRunTimeAlt());
3038       std::cout<< "reset run numbers calib: "<< calib.GetRunNumber() << "\t" << calib.GetRunNumberPed() << "\t" << calib.GetRunNumberMip() << std::endl;
3039     }
3040     double Etot = 0;
3041     int nCells  = 0;
3042     for(int j=0; j<event.GetNTiles(); j++){
3043       Caen* aTile=(Caen*)event.GetTile(j);
3044       // remove bad channels from output
3045       if (calib.GetBCCalib() && calib.GetBadChannel(aTile->GetCellID())!= 3 ){
3046         event.RemoveTile(aTile);
3047         j--;        
3048         continue;
3049       }
3050       
3051       double energy = 0;
3052       double corrHG = aTile->GetADCHigh()-calib.GetPedestalMeanH(aTile->GetCellID());
3053       double corrLG = aTile->GetADCLow()-calib.GetPedestalMeanL(aTile->GetCellID());
3054       double corrLG_HGeq = corrLG*calib.GetLGHGCorr(aTile->GetCellID()) + calib.GetLGHGCorrOff(aTile->GetCellID());
3055       if(corrHG<corrHGADCSwap){
3056         if(corrHG/calib.GetScaleHigh(aTile->GetCellID()) > minMipFrac){
3057           energy=corrHG/calib.GetScaleHigh(aTile->GetCellID());
3058         }
3059       } else {
3060         energy=corrLG/calib.GetCalcScaleLow(aTile->GetCellID());
3061       }
3062       if (debug > 1 && corrHG >= corrHGADCSwap-100 && corrHG < 4000-calib.GetPedestalMeanH(aTile->GetCellID())){
3063           std::cout << "-> Cell ID: " <<  aTile->GetCellID() << "\t HG\t" << corrHG << "\t" << corrHG/calib.GetScaleHigh(aTile->GetCellID()) << "\t LG \t" << corrLG << "\t" <<  corrLG/calib.GetCalcScaleLow(aTile->GetCellID()) << "\t"<< corrLG/calib.GetScaleLow(aTile->GetCellID()) << "\t delta: \t"<< corrHG/calib.GetScaleHigh(aTile->GetCellID())-(corrLG/calib.GetCalcScaleLow(aTile->GetCellID())) << "\tLGHG\t" << calib.GetLGHGCorr(aTile->GetCellID())<< std::endl;
3064       }
3065       // calculate trigger primitives
3066       bool localMuonTrigg   = false;
3067       bool localNoiseTrigg  = false;
3068 
3069       if (!UseLocTriggFromFile()){
3070         aTile->SetLocalTriggerPrimitive(event.CalculateLocalMuonTrigg(calib, rand, aTile->GetCellID(), localTriggerTiles, avLGHGCorr));
3071         aTile->SetLocalTriggerBit(0);
3072         localMuonTrigg   = event.InspectIfLocalMuonTrigg(aTile->GetCellID(), averageScale, factorMinTrigg, factorMaxTrigg);
3073         localNoiseTrigg  = event.InspectIfNoiseTrigg(aTile->GetCellID(), averageScale, factorMinTriggNoise);
3074         if (localMuonTrigg) aTile->SetLocalTriggerBit(1);
3075         if (localNoiseTrigg) aTile->SetLocalTriggerBit(2);
3076       } else {
3077         if (aTile->GetLocalTriggerBit() == 2)  localNoiseTrigg = true;
3078       }
3079       
3080       hspectraHGvsCellID->Fill(aTile->GetCellID(), aTile->GetADCHigh());
3081       hspectraLGvsCellID->Fill(aTile->GetCellID(), aTile->GetADCLow());
3082       hspectraHGCorrvsCellID->Fill(aTile->GetCellID(), corrHG);
3083       hspectraLGCorrvsCellID->Fill(aTile->GetCellID(), corrLG);
3084       
3085       if (localNoiseTrigg){
3086         hspectraHGCorrvsCellIDNoise->Fill(aTile->GetCellID(), corrHG);
3087         hspectraLGCorrvsCellIDNoise->Fill(aTile->GetCellID(), corrLG);
3088       
3089       }
3090       ithSpectra=hSpectra.find(aTile->GetCellID());
3091       if(ithSpectra!=hSpectra.end()){
3092         ithSpectra->second.FillExt(corrLG,corrHG,energy,corrLG_HGeq);
3093       } else {
3094         RootOutputHist->cd("IndividualCells");
3095         hSpectra[aTile->GetCellID()]=TileSpectra("Calibrate",1,aTile->GetCellID(),calib.GetTileCalib(aTile->GetCellID()),event.GetROtype(),debug);
3096         hSpectra[aTile->GetCellID()].FillExt(corrLG,corrHG,energy,corrLG_HGeq);
3097         RootOutput->cd();
3098       }
3099 
3100       if (localNoiseTrigg){
3101         ithSpectraNoise=hSpectraNoise.find(aTile->GetCellID());
3102         if(ithSpectraNoise!=hSpectraNoise.end()){
3103           ithSpectraNoise->second.FillExt(corrLG,corrHG,energy,corrLG_HGeq);
3104         } else {
3105           RootOutputHist->cd("IndividualCellsNoise");
3106           hSpectraNoise[aTile->GetCellID()]=TileSpectra("CalibrateNoise",1,aTile->GetCellID(),calib.GetTileCalib(aTile->GetCellID()),event.GetROtype(),debug);
3107           hSpectraNoise[aTile->GetCellID()].FillExt(corrLG,corrHG,energy,corrLG_HGeq);
3108           RootOutput->cd();
3109         }
3110       }
3111       
3112       if(energy!=0){ 
3113         aTile->SetE(energy);
3114         hspectraEnergyvsCellID->Fill(aTile->GetCellID(), energy);
3115         Etot=Etot+energy;
3116         nCells++;
3117       } else {
3118         event.RemoveTile(aTile);
3119         j--;
3120       }
3121     }
3122     hspectraEnergyTotvsNCells->Fill(nCells,Etot);
3123     RootOutput->cd();
3124     TdataOut->Fill();
3125   }
3126   TdataOut->Write();
3127   TsetupIn->CloneTree()->Write();
3128   
3129   if (IsCalibSaveToFile()){
3130     TString fileCalibPrint = RootOutputName;
3131     fileCalibPrint         = fileCalibPrint.ReplaceAll(".root","_calib.txt");
3132     calib.PrintCalibToFile(fileCalibPrint);
3133   }
3134   
3135   TcalibOut->Fill();
3136   TcalibOut->Write();
3137   
3138   
3139   RootOutput->Close();
3140   RootInput->Close();      
3141   
3142   RootOutputHist->cd();
3143 
3144   hspectraHGvsCellID->Write();
3145   hspectraLGvsCellID->Write();
3146   hspectraHGCorrvsCellID->Write();
3147   hspectraLGCorrvsCellID->Write();
3148   hspectraHGCorrvsCellIDNoise->Write();
3149   hspectraLGCorrvsCellIDNoise->Write();
3150   hspectraEnergyvsCellID->Write();
3151   hspectraEnergyTotvsNCells->Write();
3152   
3153   TH1D* hspectraEnergyTot = (TH1D*)hspectraEnergyTotvsNCells->ProjectionY();
3154   hspectraEnergyTot->SetDirectory(0);
3155   TH1D* hspectraNCells = (TH1D*)hspectraEnergyTotvsNCells->ProjectionX();
3156   hspectraNCells->SetDirectory(0);
3157   hspectraEnergyTot->Write("hTotEnergy");
3158   hspectraNCells->Write("hNCells");
3159 
3160   RootOutputHist->cd("IndividualCells");
3161   for(ithSpectra=hSpectra.begin(); ithSpectra!=hSpectra.end(); ++ithSpectra){
3162     ithSpectra->second.FitLGHGCorr(debug,false);
3163     ithSpectra->second.WriteExt(true);
3164   }
3165   RootOutputHist->cd("IndividualCellsNoise");
3166   for(ithSpectraNoise=hSpectraNoise.begin(); ithSpectraNoise!=hSpectraNoise.end(); ++ithSpectraNoise){
3167     ithSpectraNoise->second.WriteExt(true);
3168   }
3169   
3170   RootOutputHist->Close();
3171   //**********************************************************************
3172   //********************* Plotting ***************************************
3173   //**********************************************************************
3174   // Get run info object
3175   std::map<int,RunInfo>::iterator it=ri.find(runNr);
3176   
3177   TString outputDirPlots = GetPlotOutputDir();
3178   gSystem->Exec("mkdir -p "+outputDirPlots);
3179   
3180   //**********************************************************************
3181   // Create canvases for channel overview plotting
3182   //**********************************************************************
3183   Double_t textSizeRel = 0.035;
3184   StyleSettingsBasics("pdf");
3185   SetPlotStyle();
3186   
3187   TCanvas* canvas2DCorr = new TCanvas("canvasCorrPlots","",0,0,1450,1300);  // gives the page size
3188   DefaultCancasSettings( canvas2DCorr, 0.08, 0.13, 0.045, 0.07);
3189   canvas2DCorr->SetLogz(1);
3190   PlotSimple2D( canvas2DCorr, hspectraHGvsCellID, -10000, -10000, textSizeRel, Form("%s/HG.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
3191   PlotSimple2D( canvas2DCorr, hspectraLGvsCellID, -10000, -10000, textSizeRel, Form("%s/LG.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
3192   PlotSimple2D( canvas2DCorr, hspectraHGCorrvsCellID, -10000, -10000, textSizeRel, Form("%s/HGCorr.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
3193   PlotSimple2D( canvas2DCorr, hspectraHGCorrvsCellID, 300, -10000, textSizeRel, Form("%s/HGCorr_zoomed.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
3194   PlotSimple2D( canvas2DCorr, hspectraHGCorrvsCellIDNoise, -50, 200, -10000, textSizeRel, Form("%s/HGCorr_Noise.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true, "Local Noise triggered");
3195   PlotSimple2D( canvas2DCorr, hspectraLGCorrvsCellID, -10000, -10000, textSizeRel, Form("%s/LGCorr.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
3196   PlotSimple2D( canvas2DCorr, hspectraLGCorrvsCellID, 200, -10000, textSizeRel, Form("%s/LGCorr_zoomed.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
3197   PlotSimple2D( canvas2DCorr, hspectraLGCorrvsCellIDNoise, -50, 200, -10000, textSizeRel, Form("%s/LGCorr_Noise.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true, "Local Noise triggered");
3198   PlotSimple2D( canvas2DCorr, hspectraEnergyvsCellID, -10000, -10000, textSizeRel, Form("%s/EnergyVsCellID.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
3199   PlotSimple2D( canvas2DCorr, hspectraEnergyTotvsNCells, -10000, -10000, textSizeRel, Form("%s/EnergyTotalVsNCells.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
3200   
3201   TCanvas* canvas1DSimple = new TCanvas("canvas1DSimple","",0,0,1450,1300);  // gives the page size
3202   DefaultCancasSettings( canvas1DSimple, 0.08, 0.03, 0.03, 0.07);
3203   hspectraEnergyTot->Scale(1./evts);
3204   hspectraEnergyTot->GetYaxis()->SetTitle("counts/event");
3205   PlotSimple1D(canvas1DSimple, hspectraEnergyTot, -10000, -10000, textSizeRel, Form("%s/EnergyTot.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, Form("#LT E_{Tot} #GT = %.1f (mip/tile eq.)",hspectraEnergyTot->GetMean() ));
3206   hspectraNCells->Scale(1./evts);
3207   hspectraNCells->GetYaxis()->SetTitle("counts/event");
3208   PlotSimple1D(canvas1DSimple, hspectraNCells, -10000, -10000, textSizeRel, Form("%s/NCells.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, Form("#LT N_{Cells} #GT = %.1f",hspectraNCells->GetMean() ));
3209   
3210   if (ExtPlot > 0){
3211     //***********************************************************************************************************
3212     //********************************* 8 Panel overview plot  **************************************************
3213     //***********************************************************************************************************
3214     //*****************************************************************
3215       // Test beam geometry (beam coming from viewer)
3216       //===========================================================
3217       //||    8 (4)    ||    7 (5)   ||    6 (6)   ||    5 (7)   ||  row 0
3218       //===========================================================
3219       //||    1 (0)    ||    2 (1)   ||    3 (2)   ||    4 (3)   ||  row 1
3220       //===========================================================
3221       //    col 0     col 1       col 2     col  3
3222       // rebuild pad geom in similar way (numbering -1)
3223     //*****************************************************************
3224     TCanvas* canvas8Panel;
3225     TPad* pad8Panel[8];
3226     Double_t topRCornerX[8];
3227     Double_t topRCornerY[8];
3228     Int_t textSizePixel = 30;
3229     Double_t relSize8P[8];
3230     CreateCanvasAndPadsFor8PannelTBPlot(canvas8Panel, pad8Panel,  topRCornerX, topRCornerY, relSize8P, textSizePixel);
3231 
3232     TCanvas* canvas8PanelProf;
3233     TPad* pad8PanelProf[8];
3234     Double_t topRCornerXProf[8];
3235     Double_t topRCornerYProf[8];
3236     Double_t relSize8PProf[8];
3237     CreateCanvasAndPadsFor8PannelTBPlot(canvas8PanelProf, pad8PanelProf,  topRCornerXProf, topRCornerYProf, relSize8PProf, textSizePixel, 0.045, "Prof", false);
3238 
3239     calib.PrintGlobalInfo();  
3240     std::cout << "plotting single layers" << std::endl;
3241 
3242     for (Int_t l = 0; l < setup->GetNMaxLayer()+1; l++){    
3243       if (l%10 == 0 && l > 0 && debug > 0)
3244         std::cout << "============================== layer " <<  l << " / " << setup->GetNMaxLayer() << " layers" << std::endl;
3245       if (ExtPlot > 0){
3246         PlotNoiseAdvWithFitsFullLayer (canvas8Panel,pad8Panel, topRCornerX, topRCornerY, relSize8P, textSizePixel, 
3247                                     hSpectra, hSpectraNoise, true, -50, 100, 1.2, l, 0,
3248                                     Form("%s/NoiseTrigg_HG_Layer%02d.%s" ,outputDirPlots.Data(), l, plotSuffix.Data()), it->second);
3249         PlotNoiseAdvWithFitsFullLayer (canvas8Panel,pad8Panel, topRCornerX, topRCornerY, relSize8P, textSizePixel, 
3250                                     hSpectra, hSpectraNoise, false, -50, 100, 1.2, l, 0,
3251                                     Form("%s/NoiseTrigg_LG_Layer%02d.%s" ,outputDirPlots.Data(), l, plotSuffix.Data()), it->second);
3252         PlotSpectraFullLayer (canvas8Panel,pad8Panel, topRCornerX, topRCornerY, relSize8P, textSizePixel, 
3253                                   hSpectra, 0, -100, 4000, 1.2, l, 0,
3254                                   Form("%s/Spectra_HG_Layer%02d.%s" ,outputDirPlots.Data(), l, plotSuffix.Data()), it->second);
3255         PlotSpectraFullLayer (canvas8Panel,pad8Panel, topRCornerX, topRCornerY, relSize8P, textSizePixel, 
3256                                   hSpectra, 2, -2, 100, 1.2, l, 0,
3257                                   Form("%s/Spectra_Comb_Layer%02d.%s" ,outputDirPlots.Data(), l, plotSuffix.Data()), it->second);
3258         PlotCorrWithFitsFullLayer(canvas8PanelProf,pad8PanelProf, topRCornerXProf, topRCornerYProf, relSize8PProf, textSizePixel, 
3259                                   hSpectra, 0, -20, 800, 50., l, 0,
3260                                   Form("%s/LGHG_Corr_Layer%02d.%s" ,outputDirPlots.Data(), l, plotSuffix.Data()), it->second);
3261         PlotCorrWithFitsFullLayer(canvas8PanelProf,pad8PanelProf, topRCornerXProf, topRCornerYProf, relSize8PProf, textSizePixel, 
3262                                   hSpectra, 2, -100, 340, 300., l, 0,
3263                                   Form("%s/LGLGhgeq_Corr_Layer%02d.%s" ,outputDirPlots.Data(), l, plotSuffix.Data()), it->second);
3264       }
3265     }
3266     std::cout << "done plotting single layers" << std::endl;  
3267   }
3268   
3269   return true;
3270 }
3271 
3272 
3273 //***********************************************************************************************
3274 //*********************** Save Noise triggers only ***************************************************
3275 //***********************************************************************************************
3276 bool Analyses::SaveNoiseTriggersOnly(void){
3277   std::cout<<"Save noise triggers into separate file"<<std::endl;
3278   TcalibIn->GetEntry(0);
3279   // check whether calib should be overwritten based on external text file
3280   if (OverWriteCalib){
3281     calib.ReadCalibFromTextFile(ExternalCalibFile,debug);
3282   }
3283   
3284   int evts=TdataIn->GetEntries();
3285   for(int i=0; i<evts; i++){
3286     TdataIn->GetEntry(i);
3287     if (i%5000 == 0 && i > 0 && debug > 0) std::cout << "Reading " <<  i << " / " << evts << " events" << std::endl;
3288     for(int j=0; j<event.GetNTiles(); j++){
3289       Caen* aTile=(Caen*)event.GetTile(j);
3290       // testing for noise trigger
3291       if(aTile->GetLocalTriggerBit()!= (char)2){
3292         event.RemoveTile(aTile);
3293         j--;
3294       }
3295     }
3296     RootOutput->cd();
3297     TdataOut->Fill();
3298   }
3299   TdataOut->Write();
3300   TsetupIn->CloneTree()->Write();
3301   
3302   if (IsCalibSaveToFile()){
3303     TString fileCalibPrint = RootOutputName;
3304     fileCalibPrint         = fileCalibPrint.ReplaceAll(".root","_calib.txt");
3305     calib.PrintCalibToFile(fileCalibPrint);
3306   }
3307 
3308   TcalibOut->Fill();
3309   TcalibOut->Write();
3310   RootOutput->Close();
3311   RootInput->Close();      
3312   
3313   return true;
3314 }
3315 
3316 //***********************************************************************************************
3317 //*********************** Save Noise triggers only ***************************************************
3318 //***********************************************************************************************
3319 bool Analyses::SaveCalibToOutputOnly(void){
3320   std::cout<<"Save calib into separate file: "<< GetRootCalibOutputName().Data() <<std::endl;
3321   RootCalibOutput->cd();
3322   TcalibIn->GetEntry(0);  
3323   TsetupIn->CloneTree()->Write();
3324   
3325   // check whether calib should be overwritten based on external text file
3326   if (OverWriteCalib){
3327     calib.ReadCalibFromTextFile(ExternalCalibFile,debug);
3328   }
3329   
3330   if (IsCalibSaveToFile()){
3331     TString fileCalibPrint = GetRootCalibOutputName();
3332     fileCalibPrint         = fileCalibPrint.ReplaceAll(".root","_calib.txt");
3333     calib.PrintCalibToFile(fileCalibPrint);
3334   }
3335 
3336   TcalibOut->Fill();
3337   TcalibOut->Write();
3338   RootCalibOutput->Close();
3339   
3340   return true;
3341 }
3342 
3343 
3344 //***********************************************************************************************
3345 //*********************** Save local muon triggers only ***************************************************
3346 //***********************************************************************************************
3347 bool Analyses::SaveMuonTriggersOnly(void){
3348   std::cout<<"Save local muon triggers into separate file"<<std::endl;
3349   TcalibIn->GetEntry(0);
3350     // check whether calib should be overwritten based on external text file
3351   if (OverWriteCalib){
3352     calib.ReadCalibFromTextFile(ExternalCalibFile,debug);
3353   }
3354 
3355   int evts=TdataIn->GetEntries();
3356   for(int i=0; i<evts; i++){
3357     TdataIn->GetEntry(i);
3358     if (i%5000 == 0 && i > 0 && debug > 0) std::cout << "Reading " <<  i << " / " << evts << " events" << std::endl;
3359     for(int j=0; j<event.GetNTiles(); j++){
3360       Caen* aTile=(Caen*)event.GetTile(j);
3361       // testing for muon trigger
3362       if(aTile->GetLocalTriggerBit()!= (char)1){
3363         event.RemoveTile(aTile);
3364         j--;
3365       }
3366     }
3367     RootOutput->cd();
3368     TdataOut->Fill();
3369   }
3370   TdataOut->Write();
3371   TsetupIn->CloneTree()->Write();
3372   
3373   if (IsCalibSaveToFile()){
3374     TString fileCalibPrint = RootOutputName;
3375     fileCalibPrint         = fileCalibPrint.ReplaceAll(".root","_calib.txt");
3376     calib.PrintCalibToFile(fileCalibPrint);
3377   }
3378 
3379   TcalibOut->Fill();
3380   TcalibOut->Write();
3381   RootOutput->Close();
3382   RootInput->Close();      
3383   
3384   return true;
3385 }
3386 
3387 
3388 //***********************************************************************************************
3389 //*********************** Skim HGCROC data ******************************************************
3390 //***********************************************************************************************
3391 bool Analyses::SkimHGCROCData(void){
3392   std::cout<<"Skim HGCROC data from pure noise"<<std::endl;
3393   TcalibIn->GetEntry(0);
3394     // check whether calib should be overwritten based on external text file
3395   if (OverWriteCalib){
3396     calib.ReadCalibFromTextFile(ExternalCalibFile,debug);
3397   }
3398 
3399   int evts=TdataIn->GetEntries();
3400   if (maxEvents == -1){
3401     maxEvents = evts;
3402   } else {
3403     std::cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl;
3404     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;
3405     std::cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl;
3406   }
3407 
3408   long evtTrigg = 0;
3409   
3410   for(int i=0; i<evts && i < maxEvents; i++){
3411     TdataIn->GetEntry(i);
3412     if (i%5000 == 0 && i > 0 && debug > 0) std::cout << "Reading " <<  i << " / " << evts << " events" << std::endl;
3413     // std::cout << "Reading " <<  i << " / " << evts << " events" << std::endl;
3414     bool triggered = false;
3415     for(int j=0; j<event.GetNTiles(); j++){
3416       Hgcroc* aTile=(Hgcroc*)event.GetTile(j);
3417       // testing for any signal beyond noise
3418       // std::cout << "Cell ID:" << aTile->GetCellID() << "\t" << calib.GetPedestalMeanH(aTile->GetCellID()) << "\t" << calib.GetPedestalMeanL(aTile->GetCellID()) << "\n \tADC-wave " ;
3419       // for (int k = 0; k < (int)aTile->GetADCWaveform().size(); k++ ){
3420       //   std::cout << aTile->GetADCWaveform().at(k) << "\t" ;
3421       // }
3422       // // std::cout << "\n \tTOT-Wave ";
3423       // // for (int k = 0; k < (int)aTile->GetTOTWaveform().size(); k++ ){
3424       // //   std::cout << aTile->GetTOTWaveform().at(k) << "\t" ;
3425       // // }
3426       // std::cout << "\n \tTOA-Wave ";
3427       // for (int k = 0; k < (int)aTile->GetTOAWaveform().size(); k++ ){
3428       //   std::cout << aTile->GetTOAWaveform().at(k) << "\t" ;
3429       // }
3430       // std::cout <<"\n\t\t\t";
3431       // for (int k = 0; k < (int)aTile->GetTOAWaveform().size(); k++ )
3432       //   std::cout <<"\t";  
3433       // std::cout << " integ: "<< aTile->GetTOT() << "\t" << aTile->GetTOA() << std::endl;
3434       
3435       if (aTile->GetTOA() > 0) triggered= true;
3436       // if(aTile->GetLocalTriggerBit()!= (char)1){
3437         // event.RemoveTile(aTile);
3438         // j--;
3439       // }
3440     }
3441     
3442     if (!triggered && debug == 3){
3443       for(int j=0; j<event.GetNTiles(); j++){
3444         Hgcroc* aTile=(Hgcroc*)event.GetTile(j);
3445         // testing for any signal beyond noise
3446         std::cout << "Cell ID:" << aTile->GetCellID() << "\t" << calib.GetPedestalMeanH(aTile->GetCellID()) << "\t" << calib.GetPedestalMeanL(aTile->GetCellID()) << "\n \tADC-wave " ;
3447         for (int k = 0; k < (int)aTile->GetADCWaveform().size(); k++ ){
3448           std::cout << aTile->GetADCWaveform().at(k) << "\t" ;
3449         }
3450         std::cout << "\n \tTOT-Wave ";
3451         for (int k = 0; k < (int)aTile->GetTOTWaveform().size(); k++ ){
3452           std::cout << aTile->GetTOTWaveform().at(k) << "\t" ;
3453         }
3454         std::cout << "\n \tTOA-Wave ";
3455         for (int k = 0; k < (int)aTile->GetTOAWaveform().size(); k++ ){
3456           std::cout << aTile->GetTOAWaveform().at(k) << "\t" ;
3457         }
3458         std::cout <<"\n\t\t\t";
3459         for (int k = 0; k < (int)aTile->GetTOAWaveform().size(); k++ )
3460           std::cout <<"\t";  
3461         std::cout << " integ: "<< aTile->GetTOT() << "\t" << aTile->GetTOA() << std::endl;
3462       }
3463     }
3464     
3465     RootOutput->cd();
3466     if (triggered){
3467       evtTrigg++;
3468       TdataOut->Fill();
3469     }
3470   }
3471   TdataOut->Write();
3472   TsetupIn->CloneTree()->Write();
3473   
3474   std::cout << "Evts in: " << maxEvents << "\t skimmed: " << evtTrigg << std::endl;
3475    
3476   if (IsCalibSaveToFile()){
3477     TString fileCalibPrint = RootOutputName;
3478     fileCalibPrint         = fileCalibPrint.ReplaceAll(".root","_calib.txt");
3479     calib.PrintCalibToFile(fileCalibPrint);
3480   }
3481 
3482   TcalibOut->Fill();
3483   TcalibOut->Write();
3484   RootOutput->Close();
3485   RootInput->Close();      
3486   
3487   return true;
3488 }
3489 
3490 
3491 //***********************************************************************************************
3492 //*********************** Create output files ***************************************************
3493 //***********************************************************************************************
3494 bool Analyses::CreateOutputRootFile(void){
3495   if(Overwrite){
3496     std::cout << "overwriting exisiting output file" << std::endl;
3497     RootOutput=new TFile(RootOutputName.Data(),"RECREATE");
3498   } else{
3499     std::cout << "creating output file" << std::endl;
3500     RootOutput = new TFile(RootOutputName.Data(),"CREATE");
3501   }
3502   if(RootOutput->IsZombie()){
3503     std::cout<<"Error opening '"<<RootOutput<<"'no reachable path? Exist without force mode to overwrite?..."<<std::endl;
3504     return false;
3505   }
3506   return true;
3507 }
3508 
3509 //***********************************************************************************************
3510 //*********************** Read external bad channel map *****************************************
3511 //***********************************************************************************************
3512 std::map<int,short> Analyses::ReadExternalBadChannelMap(void){
3513   
3514   std::cout << "Reading in external mapping file" << std::endl;
3515   std::map<int,short> bcmap;
3516   
3517   std::ifstream bcmapFile;
3518   bcmapFile.open(ExternalBadChannelMap,std::ios_base::in);
3519   if (!bcmapFile) {
3520     std::cout << "ERROR: file " << ExternalBadChannelMap.Data() << " not found!" << std::endl;
3521     return bcmap;
3522   }
3523 
3524   for( TString tempLine; tempLine.ReadLine(bcmapFile, kTRUE); ) {
3525     // check if line should be considered
3526     if (tempLine.BeginsWith("%") || tempLine.BeginsWith("#")){
3527       continue;
3528     }
3529     if (debug > 1) std::cout << tempLine.Data() << std::endl;
3530 
3531     // Separate the string according to tabulators
3532     TObjArray *tempArr  = tempLine.Tokenize(" ");
3533     if(tempArr->GetEntries()<2){
3534       if (debug > 1) std::cout << "nothing to be done" << std::endl;
3535       delete tempArr;
3536       continue;
3537     } 
3538     
3539     int mod     = ((TString)((TObjString*)tempArr->At(0))->GetString()).Atoi();
3540     int layer   = ((TString)((TObjString*)tempArr->At(1))->GetString()).Atoi();
3541     int row     = ((TString)((TObjString*)tempArr->At(2))->GetString()).Atoi();
3542     int col     = ((TString)((TObjString*)tempArr->At(3))->GetString()).Atoi();
3543     short bc    = short(((TString)((TObjString*)tempArr->At(4))->GetString()).Atoi());
3544     
3545     int cellID  = setup->GetCellID( row, col, layer, mod);    
3546                 
3547     if (debug > 1) std::cout << "cellID " << cellID << "\t BC status: " << bc<< std::endl;
3548     bcmap[cellID]=bc;
3549   }
3550   std::cout << "registered " << bcmap.size() << " bad channels!" << std::endl;
3551   return bcmap;  
3552   
3553 }
[ 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!