Back to home page

EIC code displayed by LXR
 
 

    

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

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