Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-06-30 07:55:54 

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