Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-10-16 08:07:18 

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