Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-12-08 09:31:14 

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 && nADCFirst == 4 && nDiffFirstM ==1 && nDiffMaxT == 0){      // needed for summing tests
1916           aTile->PrintWaveFormDebugInfo(calib.GetPedestalMeanH(aTile->GetCellID()), calib.GetPedestalMeanL(aTile->GetCellID()), calib.GetPedestalSigL(aTile->GetCellID()));
1917           hHighestADCAbovePedVsLayer->Fill(layer,chInLayer, adc);
1918           if(ithSpectraTrigg!=hSpectraTrigg.end()){
1919             ithSpectraTrigg->second.FillExtHGCROC(adc,toa,tot,nADCFirst);
1920             ithSpectraTrigg->second.FillWaveformVsTime(aTile->GetADCWaveform(), toa, calib.GetPedestalMeanH(aTile->GetCellID()),offset);
1921           } else {
1922             RootOutputHist->cd("IndividualCellsTrigg");
1923             hSpectraTrigg[aTile->GetCellID()]=TileSpectra("signal",3,aTile->GetCellID(),calib.GetTileCalib(aTile->GetCellID()),event.GetROtype(),debug);
1924             hSpectraTrigg[aTile->GetCellID()].FillExtHGCROC(adc,toa,tot,nADCFirst);
1925             hSpectraTrigg[aTile->GetCellID()].FillWaveformVsTime(aTile->GetADCWaveform(), toa, calib.GetPedestalMeanH(aTile->GetCellID()),offset);
1926             RootOutput->cd();
1927           }
1928         }
1929       }
1930     }
1931     RootOutput->cd();
1932     TdataOut->Fill();
1933   }
1934   //RootOutput->cd();
1935   TdataOut->Write();
1936   TsetupIn->CloneTree()->Write();
1937 
1938   std::map<int,RunInfo>::iterator it=ri.find(runNr);
1939   TString outputDirPlots = GetPlotOutputDir();
1940   Double_t textSizeRel = 0.035;
1941 
1942   if ( ExtPlot > 0) {
1943     gSystem->Exec("mkdir -p "+outputDirPlots);
1944     StyleSettingsBasics("pdf");
1945     SetPlotStyle();  
1946   }
1947   
1948   if (ExtPlot > 0){
1949     
1950     TCanvas* canvas2DSigQA = new TCanvas("canvas2DSigQA","",0,0,1450,1300);  // gives the page size
1951     DefaultCancasSettings( canvas2DSigQA, 0.08, 0.13, 0.045, 0.07);
1952 
1953     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);
1954     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);
1955     PlotSimple2D( canvas2DSigQA, hSampleDiffvsCellID, -10000, setup->GetMaxCellID()+1, textSizeRel, Form("%s/SampleDiffvsCellID.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
1956     PlotSimple2D( canvas2DSigQA, hspectraTOAvsCellID, -10000, setup->GetMaxCellID()+1, textSizeRel, Form("%s/TOAvsCellID.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
1957     PlotSimple2D( canvas2DSigQA, hspectraTOTvsCellID, -10000, setup->GetMaxCellID()+1, textSizeRel, Form("%s/TOTvsCellID.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
1958     
1959     
1960     canvas2DSigQA->SetLogz(1);
1961     PlotSimple2D( canvas2DSigQA, hHighestADCAbovePedVsLayer, -10000, -10000, textSizeRel, Form("%s/MaxADCAboveNoise_vsLayer.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);    
1962     
1963     TCanvas* canvas1DSimple = new TCanvas("canvas1DSimple","",0,0,1450,1300);  // gives the page size
1964     DefaultCancasSettings( canvas1DSimple, 0.08, 0.03, 0.03, 0.07);
1965 
1966     PlotSimple1D(canvas1DSimple, hSampleTOA, -10000, (double)it->second.samples, textSizeRel, Form("%s/NSampleToA.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1);
1967     PlotSimple1D(canvas1DSimple, hSampleMaxADC, -10000, (double)it->second.samples, textSizeRel, Form("%s/NSampleMaxADC.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1);
1968     PlotSimple1D(canvas1DSimple, hSampleAboveTh, -10000, (double)it->second.samples, textSizeRel, Form("%s/NSampleAboveTh.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1);
1969     PlotSimple1D(canvas1DSimple, hSampleDiff, -10000, (double)it->second.samples, textSizeRel, Form("%s/NSampleDiff.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1);
1970     PlotSimple1D(canvas1DSimple, hSampleDiffMin, -10000, (double)it->second.samples, textSizeRel, Form("%s/NSampleDiffMin.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1);
1971 
1972     calib.PrintGlobalInfo();  
1973     
1974     std::cout << "row max: " << setup->GetNMaxRow() << "\t column max: "  << setup->GetNMaxColumn() << std::endl;
1975     
1976     if (setup->GetNMaxRow()+1 == 2 && setup->GetNMaxColumn()+1 == 4){
1977       //***********************************************************************************************************
1978       //********************************* 8 Panel overview plot  **************************************************
1979       //***********************************************************************************************************
1980       //*****************************************************************
1981         // Test beam geometry (beam coming from viewer)
1982         //===========================================================
1983         //||    8 (4)    ||    7 (5)   ||    6 (6)   ||    5 (7)   ||  row 0
1984         //===========================================================
1985         //||    1 (0)    ||    2 (1)   ||    3 (2)   ||    4 (3)   ||  row 1
1986         //===========================================================
1987         //    col 0     col 1       col 2     col  3
1988         // rebuild pad geom in similar way (numbering -1)
1989       //*****************************************************************
1990       TCanvas* canvas8Panel;
1991       TPad* pad8Panel[8];
1992       Double_t topRCornerX[8];
1993       Double_t topRCornerY[8];
1994       Int_t textSizePixel = 30;
1995       Double_t relSize8P[8];
1996       CreateCanvasAndPadsFor8PannelTBPlot(canvas8Panel, pad8Panel,  topRCornerX, topRCornerY, relSize8P, textSizePixel);
1997 
1998       TCanvas* canvas8PanelProf;
1999       TPad* pad8PanelProf[8];
2000       Double_t topRCornerXProf[8];
2001       Double_t topRCornerYProf[8];
2002       Double_t relSize8PProf[8];
2003       CreateCanvasAndPadsFor8PannelTBPlot(canvas8PanelProf, pad8PanelProf,  topRCornerXProf, topRCornerYProf, relSize8PProf, textSizePixel, 0.045, "Prof", false);
2004 
2005       std::cout << "plotting single  8M layers" << std::endl;
2006       for (Int_t l = 0; l < 5; l++){    
2007       // for (Int_t l = 0; l < setup->GetNMaxLayer()+1; l++){    
2008         for (Int_t m = 0; m < setup->GetNMaxModule()+1; m++){
2009           if (l%10 == 0 && l > 0 && debug > 0)
2010             std::cout << "============================== layer " <<  l << " / " << setup->GetNMaxLayer() << " layers" << std::endl;     
2011           PlotCorr2D8MLayer(canvas8PanelProf,pad8PanelProf, topRCornerXProf, topRCornerYProf, relSize8PProf,
2012                               textSizePixel, hSpectra, 1, -25000, (it->second.samples)*25000, 300, l, m,
2013                               Form("%s/Waveform_Mod_%02d_Layer%02d.%s" ,outputDirPlots.Data(),m,  l, plotSuffix.Data()), it->second);
2014           PlotCorr2D8MLayer(canvas8PanelProf,pad8PanelProf, topRCornerXProf, topRCornerYProf, relSize8PProf,
2015                               textSizePixel, hSpectra, 2, 0, 1024, 300, l, m,
2016                               Form("%s/TOA_ADC_Mod_%02d_Layer%02d.%s" ,outputDirPlots.Data(),m,  l, plotSuffix.Data()), it->second);
2017           PlotCorr2D8MLayer(canvas8PanelProf,pad8PanelProf, topRCornerXProf, topRCornerYProf, relSize8PProf,
2018                               textSizePixel, hSpectra, 3, 0, 1024, it->second.samples, l, m,
2019                               Form("%s/TOA_Sample_Mod_%02d_Layer%02d.%s" ,outputDirPlots.Data(),m,  l, plotSuffix.Data()), it->second);
2020           PlotCorr2D8MLayer(canvas8PanelProf,pad8PanelProf, topRCornerXProf, topRCornerYProf, relSize8PProf,
2021                             textSizePixel, hSpectraTrigg, 1, -25000, (it->second.samples)*25000, 300, l, 0,
2022                             Form("%s/WaveformSignal_Mod_%02d_Layer%02d.%s" ,outputDirPlots.Data(), m, l, plotSuffix.Data()), it->second);            
2023           PlotCorr2D8MLayer(canvas8PanelProf,pad8PanelProf, topRCornerXProf, topRCornerYProf, relSize8PProf,
2024                               textSizePixel, hSpectraTrigg, 2, 0, 1024, 300, l, m,
2025                               Form("%s/TOA_ADC_Signal_Mod_%02d_Layer%02d.%s" ,outputDirPlots.Data(),m,  l, plotSuffix.Data()), it->second);
2026           PlotCorr2D8MLayer(canvas8PanelProf,pad8PanelProf, topRCornerXProf, topRCornerYProf, relSize8PProf,
2027                               textSizePixel, hSpectraTrigg, 3, 0, 1024, it->second.samples, l, m,
2028                               Form("%s/TOA_Sample_Signal_Mod_%02d_Layer%02d.%s" ,outputDirPlots.Data(),m,  l, plotSuffix.Data()), it->second);
2029           if (ExtPlot > 1){
2030             PlotSpectra8MLayer (canvas8Panel,pad8Panel, topRCornerX, topRCornerY, relSize8P, textSizePixel, 
2031                                 hSpectra, 0, -100, 1024, 1.2, l, m,
2032                                 Form("%s/Spectra_ADC_Mod%02d_Layer%02d.%s" ,outputDirPlots.Data(), m, l, plotSuffix.Data()), it->second);
2033             PlotSpectra8MLayer (canvas8Panel,pad8Panel, topRCornerX, topRCornerY, relSize8P, textSizePixel, 
2034                                 hSpectra, 3, 0, 1024, 1.2, l, m,
2035                                 Form("%s/TOA_Mod%02d_Layer%02d.%s" ,outputDirPlots.Data(), m, l, plotSuffix.Data()), it->second);
2036             PlotSpectra8MLayer (canvas8Panel,pad8Panel, topRCornerX, topRCornerY, relSize8P, textSizePixel, 
2037                                 hSpectra, 4, 0, 4096, 1.2, l, m,
2038                                 Form("%s/TOT_Mod%02d_Layer%02d.%s" ,outputDirPlots.Data(), m, l, plotSuffix.Data()), it->second);
2039           }
2040         }
2041       }
2042       std::cout << "ending plotting single 8M layers" << std::endl;
2043     } else if (setup->GetNMaxRow()+1 == 1 && setup->GetNMaxColumn()+1 == 2){
2044       //***********************************************************************************************************
2045       //********************************* 2 Panel overview plot  **************************************************
2046       //***********************************************************************************************************
2047       //*****************************************************************
2048         // Test beam geometry (beam coming from viewer)
2049         //===========================================================
2050         //||    1 (0)    ||    2 (1)   || row 0
2051         //===========================================================
2052         //    col 0     col 1 
2053         // rebuild pad geom in similar way (numbering -1)
2054       //*****************************************************************
2055       TCanvas* canvas2Panel;
2056       TPad* pad2Panel[2];
2057       Double_t topRCornerX[2];
2058       Double_t topRCornerY[2];
2059       Int_t textSizePixel = 30;
2060       Double_t relSizeP[2];
2061       CreateCanvasAndPadsFor2PannelTBPlot(canvas2Panel, pad2Panel,  topRCornerX, topRCornerY, relSizeP, textSizePixel);
2062 
2063       TCanvas* canvas2PanelProf;
2064       TPad* pad2PanelProf[2];
2065       Double_t topRCornerXProf[2];
2066       Double_t topRCornerYProf[2];
2067       Double_t relSizePProf[2];
2068       CreateCanvasAndPadsFor2PannelTBPlot(canvas2PanelProf, pad2PanelProf,  topRCornerXProf, topRCornerYProf, relSizePProf, textSizePixel, 0.075, "Prof", false);
2069 
2070       std::cout << "plotting single  2M layers" << std::endl;
2071       for (Int_t l = 0; l < setup->GetNMaxLayer()+1; l++){    
2072         for (Int_t m = 0; m < setup->GetNMaxModule()+1; m++){
2073           if (l%10 == 0 && l > 0 && debug > 0)
2074             std::cout << "============================== layer " <<  l << " / " << setup->GetNMaxLayer() << " layers" << std::endl;     
2075           PlotCorr2D2MLayer(canvas2PanelProf,pad2PanelProf, topRCornerXProf, topRCornerYProf, relSizePProf,
2076                               textSizePixel, hSpectra, -25000, (it->second.samples)*25000, 1000, l, m,
2077                               Form("%s/Waveform_Mod_%02d_Layer%02d.%s" ,outputDirPlots.Data(),m,  l, plotSuffix.Data()), it->second);
2078           PlotCorr2D2MLayer(canvas2PanelProf,pad2PanelProf, topRCornerXProf, topRCornerYProf, relSizePProf,
2079                             textSizePixel, hSpectraTrigg, -25000, (it->second.samples)*25000, 1000, l, 0,
2080                             Form("%s/WaveformSignal_Mod_%02d_Layer%02d.%s" ,outputDirPlots.Data(), m, l, plotSuffix.Data()), it->second);            
2081           if (ExtPlot > 1){
2082             PlotSpectra2MLayer (canvas2Panel,pad2Panel, topRCornerX, topRCornerY, relSizeP, textSizePixel, 
2083                                 hSpectra, 0, -100, 1024, 1.2, l, m,
2084                                 Form("%s/Spectra_ADC_Mod%02d_Layer%02d.%s" ,outputDirPlots.Data(), m, l, plotSuffix.Data()), it->second);
2085           }          
2086         }
2087       }          
2088     } else if (setup->GetNMaxRow()+1 == 1 && setup->GetNMaxColumn()+1 == 1){
2089       //***********************************************************************************************************
2090       //********************************* Single tile overview plot  **************************************************
2091       //***********************************************************************************************************
2092       TCanvas* canvasLayer = new TCanvas("canvasLayer","",0,0,620,600);
2093       DrawCanvasSettings( canvasLayer,0.12, 0.03, 0.03, 0.1);
2094       Double_t topRCornerX = 0.95;
2095       Double_t topRCornerY = 0.95;
2096       Int_t textSizePixel = 30;
2097       Double_t relSizeP = 30./620;
2098     
2099       TCanvas* canvasLayerProf = new TCanvas("canvasLayerProf","",0,0,620,600);
2100       DrawCanvasSettings( canvasLayerProf,0.138, 0.08, 0.03, 0.1);
2101       Double_t topRCornerXProf = 0.175;
2102       Double_t topRCornerYProf = 0.95;
2103       Double_t relSizePProf = 30./620;
2104 
2105       std::cout << "plotting single tile layers" << std::endl;
2106       for (Int_t l = 0; l < setup->GetNMaxLayer()+1; l++){    
2107         for (Int_t m = 0; m < setup->GetNMaxModule()+1; m++){
2108           if (l%10 == 0 && l > 0 && debug > 0)
2109             std::cout << "============================== layer " <<  l << " / " << setup->GetNMaxLayer() << " layers" << std::endl;     
2110           PlotCorr2D1MLayer(canvasLayerProf, topRCornerXProf, topRCornerYProf, relSizePProf,
2111                               textSizePixel, hSpectra, 1, -25000, (it->second.samples)*25000, 300, l, m,
2112                               Form("%s/Waveform_Mod_%02d_Layer%02d.%s" ,outputDirPlots.Data(),m,  l, plotSuffix.Data()), it->second);
2113           PlotCorr2D1MLayer(canvasLayerProf, topRCornerXProf, topRCornerYProf, relSizePProf,
2114                             textSizePixel, hSpectra, 2, 0, 1024, 300, l, 0,
2115                             Form("%s/TOA_ADC_Mod_%02d_Layer%02d.%s" ,outputDirPlots.Data(), m, l, plotSuffix.Data()), it->second);            
2116           PlotCorr2D1MLayer(canvasLayerProf, topRCornerXProf, topRCornerYProf, relSizePProf,
2117                             textSizePixel, hSpectra, 3,0, 1024, it->second.samples, l, 0,
2118                             Form("%s/TOA_Sample_Mod_%02d_Layer%02d.%s" ,outputDirPlots.Data(), m, l, plotSuffix.Data()), it->second);            
2119           PlotCorr2D1MLayer(canvasLayerProf, topRCornerXProf, topRCornerYProf, relSizePProf,
2120                             textSizePixel, hSpectraTrigg, 1, -25000, (it->second.samples)*25000, 300, l, 0,
2121                             Form("%s/WaveformSignal_Mod_%02d_Layer%02d.%s" ,outputDirPlots.Data(), m, l, plotSuffix.Data()), it->second);            
2122           PlotCorr2D1MLayer(canvasLayerProf, topRCornerXProf, topRCornerYProf, relSizePProf,
2123                             textSizePixel, hSpectraTrigg, 2, 0, 1024, 300, l, 0,
2124                             Form("%s/TOA_ADC_Signal_Mod_%02d_Layer%02d.%s" ,outputDirPlots.Data(), m, l, plotSuffix.Data()), it->second);            
2125           PlotCorr2D1MLayer(canvasLayerProf, topRCornerXProf, topRCornerYProf, relSizePProf,
2126                             textSizePixel, hSpectraTrigg, 3,0, 1024, it->second.samples, l, 0,
2127                             Form("%s/TOA_Sample_Signal_Mod_%02d_Layer%02d.%s" ,outputDirPlots.Data(), m, l, plotSuffix.Data()), it->second);            
2128           if (ExtPlot > 1){
2129             PlotSpectra1MLayer (canvasLayer, topRCornerX, topRCornerY, relSizeP, textSizePixel, 
2130                                 hSpectra, 0, -100, 1024, 1.2, l, m,
2131                                 Form("%s/Spectra_ADC_Mod%02d_Layer%02d.%s" ,outputDirPlots.Data(), m, l, plotSuffix.Data()), it->second);
2132             PlotSpectra1MLayer (canvasLayer, topRCornerX, topRCornerY, relSizeP, textSizePixel, 
2133                                 hSpectra, 3, 0, 1024, 1.2, l, m,
2134                                 Form("%s/TOA_Mod%02d_Layer%02d.%s" ,outputDirPlots.Data(), m, l, plotSuffix.Data()), it->second);
2135             PlotSpectra1MLayer (canvasLayer, topRCornerX, topRCornerY, relSizeP, textSizePixel, 
2136                                 hSpectra, 4, 0, 4096, 1.2, l, m,
2137                                 Form("%s/TOT_Mod%02d_Layer%02d.%s" ,outputDirPlots.Data(), m, l, plotSuffix.Data()), it->second);
2138           }          
2139         }
2140       }          
2141     }
2142   }
2143   
2144   RootOutput->cd();
2145   
2146   std::cout<<"What is the value? <ped mean high>: "<<calib.GetAveragePedestalMeanHigh() << "\t good channels: " << calib.GetNumberOfChannelsWithBCflag(3) <<std::endl;
2147   if (IsCalibSaveToFile()){
2148     TString fileCalibPrint = RootOutputName;
2149     fileCalibPrint         = fileCalibPrint.ReplaceAll(".root","_calib.txt");
2150     calib.PrintCalibToFile(fileCalibPrint);
2151   }
2152   
2153   TcalibOut->Fill();
2154   TcalibOut->Write();
2155   
2156   RootOutput->Close();
2157   RootOutputHist->cd();
2158   if (typeRO == ReadOut::Type::Hgcroc){
2159     hHighestADCAbovePedVsLayer->Write();
2160     hSampleTOA->Write();
2161     hSampleMaxADC->Write();
2162     hSampleTOAVsCellID->Write();
2163     hSampleMaxADCVsCellID->Write();
2164     hspectraTOAvsCellID->Write();
2165     hspectraTOTvsCellID->Write();
2166   }
2167   RootOutputHist->cd("IndividualCells");
2168   for(ithSpectra=hSpectra.begin(); ithSpectra!=hSpectra.end(); ++ithSpectra){
2169     ithSpectra->second.WriteExt(true);
2170   }
2171   RootOutputHist->cd("IndividualCellsTrigg");
2172   if (typeRO == ReadOut::Type::Hgcroc){
2173     for(ithSpectraTrigg=hSpectraTrigg.begin(); ithSpectraTrigg!=hSpectraTrigg.end(); ++ithSpectraTrigg){
2174       ithSpectraTrigg->second.WriteExt(true);
2175     }
2176   }
2177 
2178   
2179   RootInput->Close();
2180   return true;
2181 }
2182 
2183 
2184 
2185 //***********************************************************************************************
2186 //*********************** intial scaling calculation function *********************************
2187 //***********************************************************************************************
2188 bool Analyses::GetScaling(void){
2189   std::cout<<"GetScaling"<<std::endl;
2190   
2191   std::map<int,RunInfo> ri=readRunInfosFromFile(RunListInputName.Data(),debug,0);
2192   
2193   std::map<int,TileSpectra> hSpectra;
2194   std::map<int,TileSpectra> hSpectraTrigg;
2195   std::map<int, TileSpectra>::iterator ithSpectra;
2196   std::map<int, TileSpectra>::iterator ithSpectraTrigg;
2197   
2198   ROOT::Math::MinimizerOptions::SetDefaultMinimizer("Minuit2", "Migrad");  
2199   
2200   std::cout << "Additional Output with histos being created: " << RootOutputNameHist.Data() << std::endl;
2201   if(Overwrite){
2202     std::cout << "recreating file with hists" << std::endl;
2203     RootOutputHist = new TFile(RootOutputNameHist.Data(),"RECREATE");
2204   } else{
2205     std::cout << "newly creating file with hists" << std::endl;
2206     RootOutputHist = new TFile(RootOutputNameHist.Data(),"CREATE");
2207   }
2208     
2209   //***********************************************************************************************
2210   //************************* first pass over tree to extract spectra *****************************
2211   //***********************************************************************************************
2212   // entering histoOutput file
2213   RootOutputHist->mkdir("IndividualCells");
2214   RootOutputHist->cd("IndividualCells");
2215 
2216   TcalibIn->GetEntry(0);
2217   // check whether calib should be overwritten based on external text file
2218   if (OverWriteCalib){
2219     calib.ReadCalibFromTextFile(ExternalCalibFile,debug);
2220   }
2221   
2222   int evts=TdataIn->GetEntries();
2223   int runNr = -1;
2224   if (maxEvents == -1){
2225     maxEvents = evts;
2226   } else {
2227     std::cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl;
2228     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;
2229     std::cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl;
2230   }
2231   ReadOut::Type typeRO = ReadOut::Type::Caen;
2232   int evtDeb = 5000;
2233   for(int i=0; i<evts && i < maxEvents; i++){
2234     TdataIn->GetEntry(i);
2235     if (i == 0){
2236       runNr = event.GetRunNumber();
2237       typeRO = event.GetROtype();
2238       if (typeRO != ReadOut::Type::Caen)
2239         evtDeb = 400;
2240       std::cout<< "Total number of events: " << evts << std::endl;
2241       std::cout<< "original run numbers calib: "<<calib.GetRunNumber() << "\t" << calib.GetRunNumberPed() << "\t" << calib.GetRunNumberMip() << std::endl;
2242       calib.SetRunNumberMip(runNr);
2243       calib.SetBeginRunTimeMip(event.GetBeginRunTimeAlt());
2244       std::cout<< "reset run numbers calib: "<< calib.GetRunNumber() << "\t" << calib.GetRunNumberPed() << "\t" << calib.GetRunNumberMip() << std::endl;
2245     }
2246     if (i%evtDeb == 0 && i > 0 && debug > 0) std::cout << "Reading " <<  i << " / " << evts << " events" << std::endl;
2247     if (i == 0 && debug > 2) std::cout << "decoding cell IDs" << std::endl ;
2248     for(int j=0; j<event.GetNTiles(); j++){
2249       
2250       if (typeRO == ReadOut::Type::Caen) {
2251         Caen* aTile=(Caen*)event.GetTile(j);
2252         if (i == 0 && debug > 2) std::cout << ((TString)setup->DecodeCellID(aTile->GetCellID())).Data() << std::endl;
2253         
2254         ithSpectra=hSpectra.find(aTile->GetCellID());
2255         double hgCorr = aTile->GetADCHigh()-calib.GetPedestalMeanH(aTile->GetCellID());
2256         double lgCorr = aTile->GetADCLow()-calib.GetPedestalMeanL(aTile->GetCellID());
2257             
2258         if(ithSpectra!=hSpectra.end()){
2259           ithSpectra->second.FillSpectraCAEN(lgCorr,hgCorr);
2260           if (hgCorr > 3*calib.GetPedestalSigH(aTile->GetCellID()) && lgCorr > 3*calib.GetPedestalSigL(aTile->GetCellID()) && hgCorr < 3900 )
2261             ithSpectra->second.FillCorrCAEN(lgCorr,hgCorr);
2262         } else {
2263           RootOutputHist->cd("IndividualCells");
2264           hSpectra[aTile->GetCellID()]=TileSpectra("mip1st",aTile->GetCellID(),calib.GetTileCalib(aTile->GetCellID()),event.GetROtype(), debug);
2265           hSpectra[aTile->GetCellID()].FillSpectraCAEN(lgCorr,hgCorr);;
2266           if (hgCorr > 3*calib.GetPedestalSigH(aTile->GetCellID()) && lgCorr > 3*calib.GetPedestalSigL(aTile->GetCellID() && hgCorr < 3900) )
2267             hSpectra[aTile->GetCellID()].FillCorrCAEN(lgCorr,hgCorr);
2268           RootOutput->cd();
2269         }
2270       } else if (typeRO == ReadOut::Type::Hgcroc) {
2271         Hgcroc* aTile=(Hgcroc*)event.GetTile(j);
2272         if (i == 0 && debug > 2) std::cout << ((TString)setup->DecodeCellID(aTile->GetCellID())).Data() << std::endl;
2273         
2274         ithSpectra=hSpectra.find(aTile->GetCellID());
2275         double adc = aTile->GetIntegratedADC();        
2276         double tot = aTile->GetTOT();
2277         double toa = aTile->GetTOA();
2278 
2279         if(ithSpectra!=hSpectra.end()){
2280           ithSpectra->second.FillHGCROC(adc,toa,tot);
2281         } else {
2282           RootOutputHist->cd("IndividualCells");
2283           hSpectra[aTile->GetCellID()]=TileSpectra("mip1st",aTile->GetCellID(),calib.GetTileCalib(aTile->GetCellID()),event.GetROtype(), debug);
2284           hSpectra[aTile->GetCellID()].FillHGCROC(adc,toa,tot);
2285           RootOutput->cd();
2286         }
2287         
2288       }
2289     } 
2290     RootOutput->cd();
2291   }
2292   // TdataOut->Write();
2293   TsetupIn->CloneTree()->Write();
2294   
2295   RootOutputHist->cd();
2296  
2297   //***********************************************************************************************
2298   //***** Monitoring histos for calib parameters & fits results of 1st iteration ******************
2299   //***********************************************************************************************
2300   int maxChannelPerLayer        = (setup->GetNMaxColumn()+1)*(setup->GetNMaxRow()+1);
2301   // monitoring applied pedestals
2302   TH1D* hMeanPedHGvsCellID      = new TH1D( "hMeanPedHG_vsCellID","mean Ped High Gain vs CellID ; cell ID; #mu_{noise, HG} (arb. units) ",
2303                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5);
2304   hMeanPedHGvsCellID->SetDirectory(0);
2305   TH1D* hMeanPedHG              = new TH1D( "hMeanPedHG","mean Ped High Gain ; #mu_{noise, HG} (arb. units); counts ",
2306                                             500, -0.5, 500-0.5);
2307   hMeanPedHG->SetDirectory(0);
2308   TH2D* hspectraHGMeanVsLayer   = new TH2D( "hspectraHGMeanVsLayer","Mean Ped High Gain; layer; brd channel; #mu_{Ped HG} (arb. units) ",
2309                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
2310   hspectraHGMeanVsLayer->SetDirectory(0);
2311   TH2D* hspectraHGSigmaVsLayer  = new TH2D( "hspectraHGSigmaVsLayer","Sigma Ped High Gain; layer; brd channel; #sigma_{Ped HG} (arb. units) ",
2312                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
2313   hspectraHGSigmaVsLayer->SetDirectory(0);
2314   TH1D* hMeanPedLGvsCellID      = new TH1D( "hMeanPedLG_vsCellID","mean Ped Low Gain vs CellID ; cell ID; #mu_{noise, LG} (arb. units) ",
2315                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5);
2316   hMeanPedLGvsCellID->SetDirectory(0);
2317   TH1D* hMeanPedLG             = new TH1D( "hMeanPedLG","mean Ped Low Gain ; #mu_{noise, LG} (arb. units); counts ",
2318                                             500, -0.5, 500-0.5);
2319   hMeanPedLG->SetDirectory(0);
2320   TH2D* hspectraLGMeanVsLayer   = new TH2D( "hspectraLGMeanVsLayer","Mean Ped Low Gain; layer; brd channel; #mu_{PED LG} (arb. units) ",
2321                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
2322   hspectraLGMeanVsLayer->SetDirectory(0);
2323   TH2D* hspectraLGSigmaVsLayer  = new TH2D( "hspectraLGSigmaVsLayer","Sigma Ped Low Gain; layer; brd channel; #sigma_{Ped LG} (arb. units)",
2324                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
2325   hspectraLGSigmaVsLayer->SetDirectory(0);
2326   // monitoring 1st iteration mip fits
2327   TH2D* hspectraHGMaxVsLayer1st   = new TH2D( "hspectraHGMaxVsLayer_1st","Max High Gain; layer; brd channel; Max_{HG, 1^{st}} (arb. units) ",
2328                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
2329   hspectraHGMaxVsLayer1st->SetDirectory(0);
2330   TH2D* hspectraHGFWHMVsLayer1st   = new TH2D( "hspectraHGFWHMVsLayer_1st","FWHM High Gain; layer; brd channel; FWHM_{HG, 1^{st}} (arb. units) ",
2331                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
2332   hspectraHGFWHMVsLayer1st->SetDirectory(0);
2333   TH2D* hspectraHGLMPVVsLayer1st   = new TH2D( "hspectraHGLMPVVsLayer_1st","MPV High Gain; layer; brd channel; MPV_{HG, 1^{st}} (arb. units) ",
2334                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
2335   hspectraHGLMPVVsLayer1st->SetDirectory(0);
2336   TH2D* hspectraHGLSigmaVsLayer1st = new TH2D( "hspectraHGLSigmaVsLayer_1st","Sigma Landau High Gain; layer; brd channel; #sigma_{L,HG, 1^{st}} (arb. units) ",
2337                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
2338   hspectraHGLSigmaVsLayer1st->SetDirectory(0);
2339   TH2D* hspectraHGGSigmaVsLayer1st = new TH2D( "hspectraHGGSigmaVsLayer_1st","Sigma Gauss High Gain; layer; brd channel; #sigma_{G,HG, 1^{st}} (arb. units) ",
2340                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
2341   hspectraHGGSigmaVsLayer1st->SetDirectory(0);
2342   TH2D* hLGHGCorrVsLayer = new TH2D( "hLGHGCorrVsLayer","LG-HG corr; layer; brd channel; a_{LG-HG} (arb. units) ",
2343                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
2344   hLGHGCorrVsLayer->SetDirectory(0);
2345   TH2D* hHGLGCorrVsLayer = new TH2D( "hHGLGCorrVsLayer","HG-LG corr; layer; brd channel; a_{HG-LG} (arb. units) ",
2346                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
2347   hHGLGCorrVsLayer->SetDirectory(0);
2348   TH2D* hLGHGCorrOffsetVsLayer = new TH2D( "hLGHGCorrOffsetVsLayer","LG-HG corr offset; layer; brd channel; b_{LG-HG} (arb. units) ",
2349                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
2350   hLGHGCorrOffsetVsLayer->SetDirectory(0);
2351   TH2D* hHGLGCorrOffsetVsLayer = new TH2D( "hHGLGCorrVsLayer","HG-LG corr offset; layer; brd channel; b_{HG-LG} (arb. units) ",
2352                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
2353   hHGLGCorrOffsetVsLayer->SetDirectory(0);
2354   
2355   TH1D* hMaxHG1st             = new TH1D( "hMaxHG1st","Max High Gain ;Max_{HG, 1^{st}} (arb. units) ; counts ",
2356                                             2000, -0.5, 2000-0.5);
2357   hMaxHG1st->SetDirectory(0);
2358   TH1D* hLGHGCorr             = new TH1D( "hLGHGCorr","LG-HG corr ; a_{LG-HG} (arb. units) ; counts ",
2359                                             400, -20, 20);
2360   hLGHGCorr->SetDirectory(0);
2361   TH1D* hHGLGCorr             = new TH1D( "hHGLGCorr","LG-HG corr ; a_{HG-LG} (arb. units) ; counts ",
2362                                             400, -1., 1.);
2363   hHGLGCorr->SetDirectory(0);
2364 
2365   
2366   // fit pedestal
2367   double* parameters    = new double[6];
2368   double* parErrAndRes  = new double[6];
2369   bool isGood;
2370   int currCells = 0;
2371   if ( debug > 0)
2372     std::cout << "============================== starting fitting 1st iteration" << std::endl;
2373   
2374   for(ithSpectra=hSpectra.begin(); ithSpectra!=hSpectra.end(); ++ithSpectra){
2375     if (currCells%20 == 0 && currCells > 0 && debug > 0)
2376       std::cout << "============================== cell " <<  currCells << " / " << hSpectra.size() << " cells" << std::endl;
2377     currCells++;
2378     for (int p = 0; p < 6; p++){
2379       parameters[p]   = 0;
2380       parErrAndRes[p] = 0;
2381     }
2382     isGood        = ithSpectra->second.FitMipHG(parameters, parErrAndRes, debug, yearData, false, calib.GetVov(), 1);
2383     long cellID   = ithSpectra->second.GetCellID();
2384     int layer     = setup->GetLayer(cellID);
2385     int chInLayer = setup->GetChannelInLayer(cellID);
2386     
2387     // fill cross-check histos
2388     hMeanPedHGvsCellID->SetBinContent(hMeanPedHGvsCellID->FindBin(cellID), calib.GetPedestalMeanH(cellID));
2389     hMeanPedLGvsCellID->SetBinContent(hMeanPedLGvsCellID->FindBin(cellID), calib.GetPedestalMeanL(cellID));
2390     hMeanPedHG->Fill(calib.GetPedestalMeanH(cellID));
2391     hMeanPedLG->Fill(calib.GetPedestalMeanL(cellID));
2392     
2393     int bin2D     = hspectraHGMeanVsLayer->FindBin(layer,chInLayer);
2394     hspectraHGMeanVsLayer->SetBinContent(bin2D, calib.GetPedestalMeanH(cellID));
2395     hspectraHGSigmaVsLayer->SetBinContent(bin2D, calib.GetPedestalSigH(cellID));
2396     hspectraLGMeanVsLayer->SetBinContent(bin2D, calib.GetPedestalMeanL(cellID));
2397     hspectraLGSigmaVsLayer->SetBinContent(bin2D, calib.GetPedestalSigL(cellID));
2398 
2399     if (isGood){
2400       hspectraHGMaxVsLayer1st->SetBinContent(bin2D, parameters[4]);
2401       hspectraHGFWHMVsLayer1st->SetBinContent(bin2D, parameters[5]);
2402       hspectraHGLMPVVsLayer1st->SetBinContent(bin2D, parameters[1]);
2403       hspectraHGLMPVVsLayer1st->SetBinError(bin2D, parErrAndRes[1]);
2404       hspectraHGLSigmaVsLayer1st->SetBinContent(bin2D, parameters[0]);
2405       hspectraHGLSigmaVsLayer1st->SetBinError(bin2D, parErrAndRes[0]);
2406       hspectraHGGSigmaVsLayer1st->SetBinContent(bin2D, parameters[3]);
2407       hspectraHGGSigmaVsLayer1st->SetBinError(bin2D, parErrAndRes[3]);
2408       
2409       hMaxHG1st->Fill(parameters[4]);
2410     }
2411     
2412     if (typeRO == ReadOut::Type::Caen) {
2413       isGood=ithSpectra->second.FitCorrCAEN(debug);
2414       if (ithSpectra->second.GetCorrModel(0)){
2415         hLGHGCorrVsLayer->SetBinContent(bin2D,ithSpectra->second.GetCorrModel(0)->GetParameter(1));
2416         hLGHGCorrVsLayer->SetBinError(bin2D,ithSpectra->second.GetCorrModel(0)->GetParError(1));
2417         hLGHGCorrOffsetVsLayer->SetBinContent(bin2D,ithSpectra->second.GetCorrModel(0)->GetParameter(0));
2418         hLGHGCorrOffsetVsLayer->SetBinError(bin2D,ithSpectra->second.GetCorrModel(0)->GetParError(0));
2419         hLGHGCorr->Fill(ithSpectra->second.GetCorrModel(0)->GetParameter(1));
2420       } 
2421       if (ithSpectra->second.GetCorrModel(1)){
2422         hHGLGCorrVsLayer->SetBinContent(bin2D,ithSpectra->second.GetCorrModel(1)->GetParameter(1));
2423         hHGLGCorrVsLayer->SetBinError(bin2D,ithSpectra->second.GetCorrModel(1)->GetParError(1));    
2424         hHGLGCorrOffsetVsLayer->SetBinContent(bin2D,ithSpectra->second.GetCorrModel(1)->GetParameter(0));
2425         hHGLGCorrOffsetVsLayer->SetBinError(bin2D,ithSpectra->second.GetCorrModel(1)->GetParError(0));    
2426         hHGLGCorr->Fill(ithSpectra->second.GetCorrModel(1)->GetParameter(1));
2427       }
2428     }
2429   }
2430   if ( debug > 0)
2431     std::cout << "============================== done fitting 1st iteration" << std::endl;
2432 
2433   TH1D* hHGTileSum[20];
2434   for (int c = 0; c < maxChannelPerLayer; c++ ){
2435     hHGTileSum[c] = new TH1D( Form("hHGTileSum_absChB%d",c),"av ADC surrounding cells ; ADC (arb. units); counts ",
2436                               4000, -0.5, 4000-0.5);
2437     hHGTileSum[c]->SetDirectory(0);
2438   }
2439   
2440   // setup trigger sel
2441   TRandom3* rand    = new TRandom3();
2442   Int_t localTriggerTiles = 4;
2443   double factorMinTrigg   = 0.5;
2444   double factorMaxTrigg   = 4.;
2445   if (yearData == 2023){
2446     localTriggerTiles = 6;
2447     factorMaxTrigg    = 3.;
2448   }
2449   if (typeRO == ReadOut::Type::Hgcroc){
2450     localTriggerTiles = 6;
2451     factorMinTrigg    = 0.3;
2452   }
2453   
2454   RootOutputHist->mkdir("IndividualCellsTrigg");
2455   RootOutputHist->cd("IndividualCellsTrigg");  
2456   //***********************************************************************************************
2457   //************************* first pass over tree to extract spectra *****************************
2458   //***********************************************************************************************  
2459   int actCh1st = 0;
2460   double averageScale = calib.GetAverageScaleHigh(actCh1st);
2461   double meanFWHMHG   = calib.GetAverageScaleWidthHigh();
2462   double avHGLGCorr   = calib.GetAverageHGLGCorr();
2463   double avHGLGOffCorr= calib.GetAverageHGLGCorrOff();
2464   double avLGHGCorr   = calib.GetAverageLGHGCorr();
2465   double avLGHGOffCorr= calib.GetAverageLGHGCorrOff();
2466   std::cout << "average HG mip: " << averageScale << "\t active ch: "<< actCh1st<< std::endl;
2467   for(int i=0; i<evts && i < maxEvents; i++){
2468     TdataIn->GetEntry(i);
2469     if (i%evtDeb == 0 && i > 0 && debug > 0) std::cout << "Reading " <<  i << " / " << evts << " events" << std::endl;
2470     for(int j=0; j<event.GetNTiles(); j++){
2471       if (typeRO == ReadOut::Type::Caen) {
2472         Caen* aTile=(Caen*)event.GetTile(j);
2473         if (i == 0 && debug > 2) std::cout << ((TString)setup->DecodeCellID(aTile->GetCellID())).Data() << std::endl;
2474         long currCellID = aTile->GetCellID();
2475         
2476         // read current tile
2477         ithSpectraTrigg=hSpectraTrigg.find(aTile->GetCellID());
2478         double hgCorr = aTile->GetADCHigh()-calib.GetPedestalMeanH(aTile->GetCellID());
2479         double lgCorr = aTile->GetADCLow()-calib.GetPedestalMeanL(aTile->GetCellID());
2480 
2481         aTile->SetLocalTriggerPrimitive(event.CalculateLocalMuonTrigg(calib, rand, currCellID, localTriggerTiles, avLGHGCorr));
2482         // estimate local muon trigger
2483         bool localMuonTrigg = event.InspectIfLocalMuonTrigg(currCellID, averageScale, factorMinTrigg, factorMaxTrigg);
2484         int chInLayer = setup->GetChannelInLayer(currCellID);    
2485         hHGTileSum[chInLayer]->Fill(aTile->GetLocalTriggerPrimitive());
2486         
2487         if(ithSpectraTrigg!=hSpectraTrigg.end()){
2488           ithSpectraTrigg->second.FillTrigger(aTile->GetLocalTriggerPrimitive());
2489         } else {
2490           RootOutputHist->cd("IndividualCellsTrigg");
2491           hSpectraTrigg[currCellID]=TileSpectra("mipTrigg",currCellID,calib.GetTileCalib(currCellID),event.GetROtype(),debug);
2492           hSpectraTrigg[currCellID].FillTrigger(aTile->GetLocalTriggerPrimitive());;
2493           RootOutput->cd();
2494         }
2495         // only fill tile spectra if 4 surrounding tiles on average are compatible with muon response
2496         if (localMuonTrigg){
2497           aTile->SetLocalTriggerBit(1);
2498           ithSpectraTrigg=hSpectraTrigg.find(aTile->GetCellID());
2499           ithSpectraTrigg->second.FillSpectraCAEN(lgCorr,hgCorr);
2500           if (hgCorr > 3*calib.GetPedestalSigH(currCellID) && lgCorr > 3*calib.GetPedestalSigL(currCellID) && hgCorr < 3900 )
2501             ithSpectraTrigg->second.FillCorrCAEN(lgCorr,hgCorr);
2502         }
2503       } else if (typeRO == ReadOut::Type::Hgcroc) {
2504         Hgcroc* aTile=(Hgcroc*)event.GetTile(j);
2505         if (i == 0 && debug > 2) std::cout << ((TString)setup->DecodeCellID(aTile->GetCellID())).Data() << std::endl;
2506         long currCellID = aTile->GetCellID();
2507         
2508         ithSpectraTrigg=hSpectraTrigg.find(aTile->GetCellID());
2509         // double adc = aTile->GetPedestal()+aTile->GetIntegratedADC();
2510         double adc = aTile->GetIntegratedADC();
2511         double tot = aTile->GetTOT();
2512         double toa = aTile->GetTOA();
2513 
2514         aTile->SetLocalTriggerPrimitive(event.CalculateLocalMuonTrigg(calib, rand, currCellID, localTriggerTiles, 0));
2515         // estimate local muon trigger
2516         bool localMuonTrigg = event.InspectIfLocalMuonTrigg(currCellID, averageScale, factorMinTrigg, factorMaxTrigg);
2517         int chInLayer = setup->GetChannelInLayer(currCellID);    
2518         hHGTileSum[chInLayer]->Fill(aTile->GetLocalTriggerPrimitive());
2519         
2520         if(ithSpectraTrigg!=hSpectraTrigg.end()){
2521           ithSpectraTrigg->second.FillTrigger(aTile->GetLocalTriggerPrimitive());
2522         } else {
2523           RootOutputHist->cd("IndividualCellsTrigg");
2524           hSpectraTrigg[currCellID]=TileSpectra("mipTrigg",currCellID,calib.GetTileCalib(currCellID),event.GetROtype(),debug);
2525           hSpectraTrigg[currCellID].FillTrigger(aTile->GetLocalTriggerPrimitive());;
2526           RootOutput->cd();
2527         }
2528         // only fill tile spectra if 4 surrounding tiles on average are compatible with muon response
2529         if (localMuonTrigg){
2530           aTile->SetLocalTriggerBit(1);
2531           ithSpectraTrigg=hSpectraTrigg.find(aTile->GetCellID());
2532           ithSpectraTrigg->second.FillHGCROC(adc,toa,tot);
2533           
2534         }
2535       }
2536     }
2537     RootOutput->cd();
2538     TdataOut->Fill();
2539   }
2540   TdataOut->Write();
2541   
2542   //***********************************************************************************************
2543   //***** Monitoring histos for fits results of 2nd iteration ******************
2544   //***********************************************************************************************
2545   RootOutputHist->cd();
2546   
2547   // monitoring trigger 
2548   TH2D* hmipTriggers              = new TH2D( "hmipTriggers","muon triggers; layer; brd channel; counts ",
2549                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
2550   hmipTriggers->SetDirectory(0);
2551   TH2D* hSuppresionNoise          = new TH2D( "hSuppresionNoise","S/B noise region; layer; brd channel; S/B noise region",
2552                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
2553   hSuppresionNoise->SetDirectory(0);
2554   TH2D* hSuppresionSignal         = new TH2D( "hSuppresionSignal","S/B signal region; layer; brd channel; S/B signal region",
2555                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
2556   hSuppresionSignal->SetDirectory(0);
2557 
2558   // monitoring 2nd iteration mip fits
2559   TH2D* hspectraHGMaxVsLayer2nd   = new TH2D( "hspectraHGMaxVsLayer_2nd","Max High Gain; layer; brd channel; Max_{HG, 2^{nd}} (arb. units) ",
2560                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
2561   hspectraHGMaxVsLayer2nd->SetDirectory(0);
2562   TH2D* hspectraHGFWHMVsLayer2nd   = new TH2D( "hspectraHGFWHMVsLayer_2nd","FWHM High Gain; layer; brd channel; FWHM_{HG, 2^{nd}} (arb. units) ",
2563                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
2564   hspectraHGFWHMVsLayer2nd->SetDirectory(0);
2565   TH2D* hspectraHGLMPVVsLayer2nd   = new TH2D( "hspectraHGLMPVVsLayer_2nd","MPV High Gain; layer; brd channel; MPV_{HG, 2^{nd}} (arb. units) ",
2566                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
2567   hspectraHGLMPVVsLayer2nd->SetDirectory(0);
2568   TH2D* hspectraHGLSigmaVsLayer2nd = new TH2D( "hspectraHGLSigmaVsLayer_2nd","Sigma Landau High Gain; layer; brd channel; #sigma_{L,HG, 2^{nd}} (arb. units) ",
2569                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
2570   hspectraHGLSigmaVsLayer2nd->SetDirectory(0);
2571   TH2D* hspectraHGGSigmaVsLayer2nd = new TH2D( "hspectraHGGSigmaVsLayer_2nd","Sigma Gauss High Gain; layer; brd channel; #sigma_{G,HG, 2^{nd}} (arb. units) ",
2572                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
2573   hspectraHGGSigmaVsLayer2nd->SetDirectory(0);
2574   TH2D* hspectraLGMaxVsLayer2nd   = new TH2D( "hspectraLGMaxVsLayer_2nd","Max High Gain; layer; brd channel; Max_{LG, 2^{nd}} (arb. units) ",
2575                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
2576   hspectraLGMaxVsLayer2nd->SetDirectory(0);
2577   TH2D* hspectraLGFWHMVsLayer2nd   = new TH2D( "hspectraLGFWHMVsLayer_2nd","FWHM High Gain; layer; brd channel; FWHM_{LG, 2^{nd}} (arb. units) ",
2578                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
2579   hspectraLGFWHMVsLayer2nd->SetDirectory(0);
2580   TH2D* hspectraLGLMPVVsLayer2nd   = new TH2D( "hspectraLGLMPVVsLayer_2nd","MPV High Gain; layer; brd channel; MPV_{LG, 2^{nd}} (arb. units) ",
2581                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
2582   hspectraLGLMPVVsLayer2nd->SetDirectory(0);
2583   TH2D* hspectraLGLSigmaVsLayer2nd = new TH2D( "hspectraLGLSigmaVsLayer_2nd","Sigma Landau High Gain; layer; brd channel; #sigma_{L,LG, 2^{nd}} (arb. units) ",
2584                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
2585   hspectraLGLSigmaVsLayer2nd->SetDirectory(0);
2586   TH2D* hspectraLGGSigmaVsLayer2nd = new TH2D( "hspectraLGGSigmaVsLayer_2nd","Sigma Gauss High Gain; layer; brd channel; #sigma_{G,LG, 2^{nd}} (arb. units) ",
2587                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
2588   hspectraLGGSigmaVsLayer2nd->SetDirectory(0);
2589 
2590   TH1D* hMaxHG2nd             = new TH1D( "hMaxHG2nd","Max High Gain ;Max_{HG, 2^{nd}} (arb. units) ; counts ",
2591                                             2000, -0.5, 2000-0.5);
2592   hMaxHG2nd->SetDirectory(0);
2593   TH1D* hMaxLG2nd             = new TH1D( "hMaxLG2nd","Max High Gain ;Max_{LG, 2^{nd}} (arb. units) ; counts ",
2594                                             400, -0.5, 400-0.5);
2595   hMaxLG2nd->SetDirectory(0);
2596 
2597 
2598   currCells = 0;
2599   if ( debug > 0)
2600     std::cout << "============================== starting fitting 2nd iteration" << std::endl;
2601   for(ithSpectraTrigg=hSpectraTrigg.begin(); ithSpectraTrigg!=hSpectraTrigg.end(); ++ithSpectraTrigg){
2602     if (currCells%20 == 0 && currCells > 0 && debug > 0)
2603       std::cout << "============================== cell " <<  currCells << " / " << hSpectraTrigg.size() << " cells" << std::endl;
2604     currCells++;
2605     for (int p = 0; p < 6; p++){
2606       parameters[p]   = 0;
2607       parErrAndRes[p] = 0;
2608     }
2609     isGood=ithSpectraTrigg->second.FitMipHG(parameters, parErrAndRes, debug, yearData, false, calib.GetVov(), averageScale);
2610     
2611     long cellID     = ithSpectraTrigg->second.GetCellID();
2612     int layer       = setup->GetLayer(cellID);
2613     int chInLayer   = setup->GetChannelInLayer(cellID);    
2614     int bin2D       = hspectraHGMeanVsLayer->FindBin(layer,chInLayer);
2615 
2616     double pedSigHG = 0;
2617     double maxBin   = 0;
2618     if (typeRO == ReadOut::Type::Caen){
2619       pedSigHG = calib.GetPedestalSigH(cellID);
2620       maxBin   = 3800;
2621     } else {
2622       pedSigHG = 20;
2623       maxBin   = 1024;
2624     }
2625     Int_t binNLow   = ithSpectraTrigg->second.GetHG()->FindBin(-1*pedSigHG);
2626     Int_t binNHigh  = ithSpectraTrigg->second.GetHG()->FindBin(3*pedSigHG);
2627     Int_t binSHigh  = ithSpectraTrigg->second.GetHG()->FindBin(maxBin);
2628     
2629     double S_NoiseR = ithSpectraTrigg->second.GetHG()->Integral(binNLow, binNHigh);
2630     double S_SigR   = ithSpectraTrigg->second.GetHG()->Integral(binNHigh, binSHigh);
2631     
2632     ithSpectra      = hSpectra.find(cellID);
2633     double B_NoiseR = ithSpectra->second.GetHG()->Integral(binNLow , binNHigh);
2634     double B_SigR   = ithSpectra->second.GetHG()->Integral(binNHigh, binSHigh);
2635     
2636     double SB_NoiseR  = (B_NoiseR != 0.) ? S_NoiseR/B_NoiseR : 0;
2637     double SB_SigR    = (B_SigR != 0.) ? S_SigR/B_SigR : 0;
2638     
2639     hmipTriggers->SetBinContent(bin2D, ithSpectraTrigg->second.GetHG()->GetEntries());
2640     hSuppresionNoise->SetBinContent(bin2D, SB_NoiseR);
2641     hSuppresionSignal->SetBinContent(bin2D, SB_SigR);
2642     if (isGood){
2643       hspectraHGMaxVsLayer2nd->SetBinContent(bin2D, parameters[4]);
2644       hspectraHGFWHMVsLayer2nd->SetBinContent(bin2D, parameters[5]);
2645       hspectraHGLMPVVsLayer2nd->SetBinContent(bin2D, parameters[1]);
2646       hspectraHGLMPVVsLayer2nd->SetBinError(bin2D, parErrAndRes[1]);
2647       hspectraHGLSigmaVsLayer2nd->SetBinContent(bin2D, parameters[0]);
2648       hspectraHGLSigmaVsLayer2nd->SetBinError(bin2D, parErrAndRes[0]);
2649       hspectraHGGSigmaVsLayer2nd->SetBinContent(bin2D, parameters[3]);
2650       hspectraHGGSigmaVsLayer2nd->SetBinError(bin2D, parErrAndRes[3]);
2651       hMaxHG2nd->Fill(parameters[4]);
2652     }
2653     
2654     if (typeRO == ReadOut::Type::Caen) {
2655       for (int p = 0; p < 6; p++){
2656         parameters[p]   = 0;
2657         parErrAndRes[p] = 0;
2658       }
2659       isGood=ithSpectraTrigg->second.FitMipLG(parameters, parErrAndRes, debug, yearData, false, 1);
2660       if (isGood){
2661         hspectraLGMaxVsLayer2nd->SetBinContent(bin2D, parameters[4]);
2662         hspectraLGFWHMVsLayer2nd->SetBinContent(bin2D, parameters[5]);
2663         hspectraLGLMPVVsLayer2nd->SetBinContent(bin2D, parameters[1]);
2664         hspectraLGLMPVVsLayer2nd->SetBinError(bin2D, parErrAndRes[1]);
2665         hspectraLGLSigmaVsLayer2nd->SetBinContent(bin2D, parameters[0]);
2666         hspectraLGLSigmaVsLayer2nd->SetBinError(bin2D, parErrAndRes[0]);
2667         hspectraLGGSigmaVsLayer2nd->SetBinContent(bin2D, parameters[3]);
2668         hspectraLGGSigmaVsLayer2nd->SetBinError(bin2D, parErrAndRes[3]);
2669         hMaxLG2nd->Fill(parameters[4]);
2670       }
2671     }
2672   }
2673   if ( debug > 0)
2674     std::cout << "============================== done fitting 2nd iteration" << std::endl;
2675   int actCh2nd = 0;
2676   double averageScaleUpdated    = calib.GetAverageScaleHigh(actCh2nd);
2677   double meanFWHMHGUpdated      = calib.GetAverageScaleWidthHigh();
2678   double averageScaleLowUpdated = 0.;
2679   double meanFWHMLGUpdated      = 0.;
2680   if (typeRO == ReadOut::Type::Caen){
2681     averageScaleLowUpdated = calib.GetAverageScaleLow();
2682     meanFWHMLGUpdated      = calib.GetAverageScaleWidthLow();
2683   }
2684   std::cout << "average 1st iteration HG mip: " << averageScale << "\t act. channels: " <<   actCh1st << std::endl;
2685   std::cout << "average 2nd iteration  HG mip: " << averageScaleUpdated << "\t LG mip: " << averageScaleLowUpdated << "\t act. channels: " <<   actCh2nd << std::endl;
2686   
2687   RootOutput->cd();
2688   
2689   if (IsCalibSaveToFile()){
2690     TString fileCalibPrint = RootOutputName;
2691     fileCalibPrint         = fileCalibPrint.ReplaceAll(".root","_calib.txt");
2692     calib.PrintCalibToFile(fileCalibPrint);
2693   }
2694 
2695   TcalibOut->Fill();
2696   TcalibOut->Write();
2697   RootOutput->Close();
2698 
2699 
2700   RootOutputHist->cd("IndividualCells");
2701     for(ithSpectra=hSpectra.begin(); ithSpectra!=hSpectra.end(); ++ithSpectra){
2702       ithSpectra->second.Write(true);
2703     }
2704   RootOutputHist->cd("IndividualCellsTrigg");
2705     for(ithSpectra=hSpectraTrigg.begin(); ithSpectra!=hSpectraTrigg.end(); ++ithSpectra){
2706       ithSpectra->second.Write(true);
2707     }
2708   RootOutputHist->cd();
2709     hMeanPedHGvsCellID->Write();
2710     hMeanPedHG->Write();
2711     
2712     hspectraHGMeanVsLayer->Write();
2713     hspectraHGSigmaVsLayer->Write();
2714     hspectraHGMaxVsLayer1st->Write();
2715     hspectraHGFWHMVsLayer1st->Write();
2716     hspectraHGLMPVVsLayer1st->Write();
2717     hspectraHGLSigmaVsLayer1st->Write();
2718     hspectraHGGSigmaVsLayer1st->Write();
2719     hMaxHG1st->Write();
2720     
2721     
2722     hspectraHGMaxVsLayer2nd->Write();
2723     hspectraHGFWHMVsLayer2nd->Write();
2724     hspectraHGLMPVVsLayer2nd->Write();
2725     hspectraHGLSigmaVsLayer2nd->Write();
2726     hspectraHGGSigmaVsLayer2nd->Write();
2727     hMaxHG2nd->Write();
2728     
2729     if (typeRO == ReadOut::Type::Caen) {
2730       hMeanPedLGvsCellID->Write();
2731       hMeanPedLG->Write();
2732       hspectraLGMeanVsLayer->Write();
2733       hspectraLGSigmaVsLayer->Write();
2734       hLGHGCorrVsLayer->Write();
2735       hLGHGCorrOffsetVsLayer->Write();
2736       hHGLGCorrVsLayer->Write();
2737       hHGLGCorrOffsetVsLayer->Write();
2738       hLGHGCorr->Write();
2739       hHGLGCorr->Write();
2740       hspectraLGMaxVsLayer2nd->Write();
2741       hspectraLGFWHMVsLayer2nd->Write();
2742       hspectraLGLMPVVsLayer2nd->Write();
2743       hspectraLGLSigmaVsLayer2nd->Write();
2744       hspectraLGGSigmaVsLayer2nd->Write();
2745       hMaxLG2nd->Write();
2746     }
2747     for (int c = 0; c< maxChannelPerLayer; c++)
2748       hHGTileSum[c]->Write();
2749     hmipTriggers->Write();
2750     hSuppresionNoise->Write();
2751     hSuppresionSignal->Write();
2752   // fill calib tree & write it
2753   // close open root files
2754   RootOutputHist->Write();
2755   RootOutputHist->Close();
2756 
2757   RootInput->Close();
2758 
2759   // Get run info object
2760   std::map<int,RunInfo>::iterator it=ri.find(runNr);
2761   
2762   // create directory for plot output
2763   TString outputDirPlots = GetPlotOutputDir();
2764   gSystem->Exec("mkdir -p "+outputDirPlots);
2765   
2766   //**********************************************************************
2767   // Create canvases for channel overview plotting
2768   //**********************************************************************
2769   Double_t textSizeRel = 0.035;
2770   StyleSettingsBasics("pdf");
2771   SetPlotStyle();
2772 
2773   TCanvas* canvas2DCorr = new TCanvas("canvasCorrPlots","",0,0,1450,1300);  // gives the page size
2774   DefaultCancasSettings( canvas2DCorr, 0.08, 0.13, 0.045, 0.07);
2775 
2776   canvas2DCorr->SetLogz(0);
2777   PlotSimple2D( canvas2DCorr, hspectraHGMeanVsLayer, -10000, -10000, textSizeRel, Form("%s/HG_NoiseMean.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
2778   PlotSimple2D( canvas2DCorr, hspectraHGSigmaVsLayer,-10000, -10000, textSizeRel, Form("%s/HG_NoiseSigma.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1,  kFALSE, "colz", true);
2779   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));
2780   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));
2781   PlotSimple2D( canvas2DCorr, hspectraHGLMPVVsLayer1st, -10000, -10000, textSizeRel, Form("%s/HG_LandMPVMip_1st.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
2782   PlotSimple2D( canvas2DCorr, hspectraHGLSigmaVsLayer1st, -10000, -10000, textSizeRel, Form("%s/HG_LandSigMip_1st.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
2783   PlotSimple2D( canvas2DCorr, hspectraHGGSigmaVsLayer1st, -10000, -10000, textSizeRel, Form("%s/HG_GaussSigMip_1st.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
2784   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));
2785   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));
2786   PlotSimple2D( canvas2DCorr, hspectraHGLMPVVsLayer2nd, -10000, -10000, textSizeRel, Form("%s/HG_LandMPVMip_2nd.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
2787   PlotSimple2D( canvas2DCorr, hspectraHGLSigmaVsLayer2nd, -10000, -10000, textSizeRel, Form("%s/HG_LandSigMip_2nd.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
2788   PlotSimple2D( canvas2DCorr, hspectraHGGSigmaVsLayer2nd, -10000, -10000, textSizeRel, Form("%s/HG_GaussSigMip_2nd.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
2789   canvas2DCorr->SetLogz(1);
2790   TString drawOpt = "colz";
2791   if (yearData == 2023)
2792     drawOpt = "colz,text";
2793   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));
2794   PlotSimple2D( canvas2DCorr, hSuppresionNoise, -10000, -10000, textSizeRel, Form("%s/SuppressionNoise.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, drawOpt, true);
2795   PlotSimple2D( canvas2DCorr, hSuppresionSignal, -10000, -10000, textSizeRel, Form("%s/SuppressionSignal.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, drawOpt, true);
2796   
2797   canvas2DCorr->SetLogz(0);
2798   
2799   if (typeRO == ReadOut::Type::Caen) {
2800     PlotSimple2D( canvas2DCorr, hspectraLGMeanVsLayer, -10000, -10000, textSizeRel, Form("%s/LG_NoiseMean.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
2801     PlotSimple2D( canvas2DCorr, hspectraLGSigmaVsLayer, -10000, -10000, textSizeRel, Form("%s/LG_NoiseSigma.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
2802     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));
2803     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));
2804     PlotSimple2D( canvas2DCorr, hspectraLGLMPVVsLayer2nd, -10000, -10000, textSizeRel, Form("%s/LG_LandMPVMip_2nd.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
2805     PlotSimple2D( canvas2DCorr, hspectraLGLSigmaVsLayer2nd, -10000, -10000, textSizeRel, Form("%s/LG_LandSigMip_2nd.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
2806     PlotSimple2D( canvas2DCorr, hspectraLGGSigmaVsLayer2nd, -10000, -10000, textSizeRel, Form("%s/LG_GaussSigMip_2nd.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
2807 
2808     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));
2809     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));
2810     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));
2811     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));
2812   }
2813   if (ExtPlot > 0){
2814     //***********************************************************************************************************
2815     //********************************* 8 Panel overview plot  **************************************************
2816     //***********************************************************************************************************
2817     //*****************************************************************
2818       // Test beam geometry (beam coming from viewer)
2819       //===========================================================
2820       //||    8 (4)    ||    7 (5)   ||    6 (6)   ||    5 (7)   ||  row 0
2821       //===========================================================
2822       //||    1 (0)    ||    2 (1)   ||    3 (2)   ||    4 (3)   ||  row 1
2823       //===========================================================
2824       //    col 0     col 1       col 2     col  3
2825       // rebuild pad geom in similar way (numbering -1)
2826     //*****************************************************************
2827     TCanvas* canvas8Panel;
2828     TPad* pad8Panel[8];
2829     Double_t topRCornerX[8];
2830     Double_t topRCornerY[8];
2831     Int_t textSizePixel = 30;
2832     Double_t relSize8P[8];
2833     CreateCanvasAndPadsFor8PannelTBPlot(canvas8Panel, pad8Panel,  topRCornerX, topRCornerY, relSize8P, textSizePixel);
2834 
2835     TCanvas* canvas8PanelProf;
2836     TPad* pad8PanelProf[8];
2837     Double_t topRCornerXProf[8];
2838     Double_t topRCornerYProf[8];
2839     Double_t relSize8PProf[8];
2840     CreateCanvasAndPadsFor8PannelTBPlot(canvas8PanelProf, pad8PanelProf,  topRCornerXProf, topRCornerYProf, relSize8PProf, textSizePixel, 0.045, "Prof", false);
2841 
2842     calib.PrintGlobalInfo();  
2843     Double_t maxHG = ReturnMipPlotRangeDepVov(calib.GetVov(),true, typeRO);
2844     Double_t maxLG = ReturnMipPlotRangeDepVov(calib.GetVov(),false, typeRO);
2845     std::cout << "plotting single layers" << std::endl;
2846 
2847     for (Int_t l = 0; l < setup->GetNMaxLayer()+1; l++){
2848       for (Int_t m = 0; m < setup->GetNMaxModule()+1; m++){    
2849         if (l%10 == 0 && l > 0 && debug > 0)
2850           std::cout << "============================== layer " <<  l << " / " << setup->GetNMaxLayer() << " layers" << std::endl;
2851         if (ExtPlot > 0){
2852           PlotMipWithFits8MLayer (canvas8Panel,pad8Panel, topRCornerX, topRCornerY, relSize8P, textSizePixel, 
2853                                     hSpectra, hSpectraTrigg, setup, true, -100, maxHG, 1.2, l, m,
2854                                     Form("%s/MIP_HG_Mod%02d_Layer%02d.%s" ,outputDirPlots.Data(), m, l, plotSuffix.Data()), it->second);
2855           Double_t maxTriggPPlot = maxHG*2;
2856           if (typeRO != ReadOut::Type::Caen)
2857             maxTriggPPlot = 500;
2858 
2859           PlotTriggerPrimWithFits8MLayer (canvas8Panel,pad8Panel, topRCornerX, topRCornerY, relSize8P, textSizePixel, 
2860                                             hSpectraTrigg, setup, averageScale, factorMinTrigg, factorMaxTrigg,
2861                                             0, maxTriggPPlot, 1.2, l, m, Form("%s/TriggPrimitive_Mod%02d_Layer%02d.%s" ,outputDirPlots.Data(), m, l, plotSuffix.Data()), it->second);
2862           if (typeRO == ReadOut::Type::Caen) {
2863             PlotCorrWithFits8MLayer(canvas8PanelProf,pad8PanelProf, topRCornerXProf, topRCornerYProf, relSize8PProf, textSizePixel, 
2864                                     hSpectra, 0, -20, 800, 3900, l, m,
2865                                     Form("%s/LGHG_Corr_Mod%02d_Layer%02d.%s" ,outputDirPlots.Data(), m, l, plotSuffix.Data()), it->second);
2866           }
2867         }
2868         if (ExtPlot > 1 && typeRO == ReadOut::Type::Caen) {
2869           PlotMipWithFits8MLayer (canvas8Panel,pad8Panel, topRCornerX, topRCornerY, relSize8P, textSizePixel, 
2870                                     hSpectra, hSpectraTrigg, setup, false, -30, maxLG, 1.2, l, m,
2871                                     Form("%s/MIP_LG_Mod%02d_Layer%02d.%s" ,outputDirPlots.Data(), m, l, plotSuffix.Data()), it->second);
2872           PlotCorrWithFits8MLayer(canvas8PanelProf,pad8PanelProf, topRCornerXProf, topRCornerYProf, relSize8PProf, textSizePixel, 
2873                                     hSpectra, 1, -100, 4000, 340, l, m,
2874                                     Form("%s/HGLG_Corr_Mod%02d_Layer%02d.%s" ,outputDirPlots.Data(), m, l, plotSuffix.Data()), it->second);
2875         }
2876       }
2877     }
2878     std::cout << "done plotting single layers" << std::endl;  
2879   }
2880   return true;
2881 }
2882 
2883 //***********************************************************************************************
2884 //*********************** improved scaling calculation function *********************************
2885 //***********************************************************************************************
2886 bool Analyses::GetImprovedScaling(void){
2887   std::cout<<"GetImprovedScaling"<<std::endl;
2888   
2889   std::map<int,RunInfo> ri=readRunInfosFromFile(RunListInputName.Data(),debug,0);
2890   
2891   std::map<int,TileSpectra> hSpectra;
2892   std::map<int,TileSpectra> hSpectraTrigg;
2893   std::map<int, TileSpectra>::iterator ithSpectra;
2894   std::map<int, TileSpectra>::iterator ithSpectraTrigg;
2895   
2896   std::cout << "Additional Output with histos being created: " << RootOutputNameHist.Data() << std::endl;
2897   if(Overwrite){
2898     std::cout << "recreating file with hists" << std::endl;
2899     RootOutputHist = new TFile(RootOutputNameHist.Data(),"RECREATE");
2900   } else{
2901     std::cout << "newly creating file with hists" << std::endl;
2902     RootOutputHist = new TFile(RootOutputNameHist.Data(),"CREATE");
2903   }
2904     
2905   ROOT::Math::MinimizerOptions::SetDefaultMinimizer("Minuit2", "Migrad");
2906   // setup trigger sel
2907   double factorMinTrigg   = 0.8;
2908   double factorMaxTrigg   = 2.5;
2909   if (yearData == 2023){
2910     factorMinTrigg    = 0.9;
2911     factorMaxTrigg    = 2.;
2912   }
2913   
2914   RootOutputHist->mkdir("IndividualCells");
2915   RootOutputHist->mkdir("IndividualCellsTrigg");
2916   RootOutputHist->cd("IndividualCellsTrigg");  
2917   //***********************************************************************************************
2918   //************************* first pass over tree to extract spectra *****************************
2919   //***********************************************************************************************  
2920   TcalibIn->GetEntry(0);
2921   // check whether calib should be overwritten based on external text file
2922   if (OverWriteCalib){
2923     calib.ReadCalibFromTextFile(ExternalCalibFile,debug);
2924   }
2925   
2926   int evts=TdataIn->GetEntries();
2927   int runNr = -1;
2928   int actChI  = 0;
2929   ReadOut::Type typeRO = ReadOut::Type::Caen;
2930   int evtDeb = 5000;
2931   
2932   if (maxEvents == -1){
2933     maxEvents = evts;
2934   } else {
2935     std::cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl;
2936     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;
2937     std::cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl;
2938   }
2939   
2940   double averageScale     = calib.GetAverageScaleHigh(actChI);
2941   double averageScaleLow  = calib.GetAverageScaleLow();
2942   std::cout << "average HG mip: " << averageScale << " LG mip: "<< averageScaleLow << "\t act. ch: "<< actChI << std::endl;
2943   
2944   for(int i=0; i<maxEvents; i++){
2945     TdataIn->GetEntry(i);    
2946     if (i == 0){
2947       runNr = event.GetRunNumber();
2948       typeRO = event.GetROtype();
2949       if (typeRO != ReadOut::Type::Caen){
2950         evtDeb = 400;
2951         factorMinTrigg    = 0.5;
2952         std::cout << "reseting lower trigger factor limit to: " << factorMinTrigg << std::endl;
2953       }
2954       std::cout<< "Total number of events: " << evts << std::endl;
2955       std::cout<< "original run numbers calib: "<<calib.GetRunNumber() << "\t" << calib.GetRunNumberPed() << "\t" << calib.GetRunNumberMip() << std::endl;
2956     }
2957 
2958     TdataIn->GetEntry(i);
2959     if (i%evtDeb == 0 && i > 0 && debug > 0) std::cout << "Reading " <<  i << " / " << maxEvents << " events" << std::endl;
2960     for(int j=0; j<event.GetNTiles(); j++){
2961       // CAEN treatment
2962       if (typeRO == ReadOut::Type::Caen) {
2963         Caen* aTile=(Caen*)event.GetTile(j);
2964         if (i == 0 && debug > 2) std::cout << ((TString)setup->DecodeCellID(aTile->GetCellID())).Data() << std::endl;
2965         long currCellID = aTile->GetCellID();
2966         
2967         // read current tile
2968         ithSpectraTrigg=hSpectraTrigg.find(currCellID);
2969         double hgCorr = aTile->GetADCHigh()-calib.GetPedestalMeanH(currCellID);
2970         double lgCorr = aTile->GetADCLow()-calib.GetPedestalMeanL(currCellID);
2971 
2972         // estimate local muon trigger
2973         bool localMuonTrigg = event.InspectIfLocalMuonTrigg(currCellID, averageScale, factorMinTrigg, factorMaxTrigg);
2974         
2975         if(ithSpectraTrigg!=hSpectraTrigg.end()){
2976           ithSpectraTrigg->second.FillTrigger(aTile->GetLocalTriggerPrimitive());
2977         } else {
2978           RootOutputHist->cd("IndividualCellsTrigg");
2979           hSpectraTrigg[currCellID]=TileSpectra("mipTrigg",currCellID,calib.GetTileCalib(currCellID),event.GetROtype(),debug);
2980           hSpectraTrigg[currCellID].FillTrigger(aTile->GetLocalTriggerPrimitive());;
2981           RootOutput->cd();
2982         }
2983         
2984         ithSpectra=hSpectra.find(currCellID);
2985         if (ithSpectra!=hSpectra.end()){
2986           ithSpectra->second.FillSpectraCAEN(lgCorr,hgCorr);
2987           if (hgCorr > 3*calib.GetPedestalSigH(currCellID) && lgCorr > 3*calib.GetPedestalSigL(currCellID) && hgCorr < 3900 )
2988             ithSpectra->second.FillCorrCAEN(lgCorr,hgCorr);
2989         } else {
2990           RootOutputHist->cd("IndividualCells");
2991           hSpectra[currCellID]=TileSpectra("mip1st",currCellID,calib.GetTileCalib(currCellID),event.GetROtype(),debug);
2992           hSpectra[currCellID].FillSpectraCAEN(lgCorr,hgCorr);;
2993           if (hgCorr > 3*calib.GetPedestalSigH(currCellID) && lgCorr > 3*calib.GetPedestalSigL(currCellID && hgCorr < 3900) )
2994             hSpectra[currCellID].FillCorrCAEN(lgCorr,hgCorr);;
2995 
2996           RootOutput->cd();
2997         }
2998         
2999       
3000         // only fill tile spectra if 4 surrounding tiles on average are compatible with muon response
3001         if (localMuonTrigg){
3002           aTile->SetLocalTriggerBit(1);
3003           ithSpectraTrigg=hSpectraTrigg.find(aTile->GetCellID());
3004           ithSpectraTrigg->second.FillSpectraCAEN(lgCorr,hgCorr);
3005           if (hgCorr > 3*calib.GetPedestalSigH(currCellID) && lgCorr > 3*calib.GetPedestalSigL(currCellID) && hgCorr < 3900 )
3006             ithSpectraTrigg->second.FillCorrCAEN(lgCorr,hgCorr);
3007         }
3008       // HGCROC treatment
3009       } else if (typeRO == ReadOut::Type::Hgcroc) {
3010         Hgcroc* aTile=(Hgcroc*)event.GetTile(j);
3011         if (i == 0 && debug > 2) std::cout << ((TString)setup->DecodeCellID(aTile->GetCellID())).Data() << std::endl;
3012         long currCellID = aTile->GetCellID();
3013 
3014         ithSpectraTrigg=hSpectraTrigg.find(currCellID);
3015         // double adc = aTile->GetPedestal()+aTile->GetIntegratedADC();
3016         double adc = aTile->GetIntegratedADC();
3017         double tot = aTile->GetTOT();
3018         double toa = aTile->GetTOA();
3019         
3020         // estimate local muon trigger
3021         bool localMuonTrigg = event.InspectIfLocalMuonTrigg(currCellID, averageScale, factorMinTrigg, factorMaxTrigg);
3022         
3023         if(ithSpectraTrigg!=hSpectraTrigg.end()){
3024           ithSpectraTrigg->second.FillTrigger(aTile->GetLocalTriggerPrimitive());
3025         } else {
3026           RootOutputHist->cd("IndividualCellsTrigg");
3027           hSpectraTrigg[currCellID]=TileSpectra("mipTrigg",currCellID,calib.GetTileCalib(currCellID),event.GetROtype(),debug);
3028           hSpectraTrigg[currCellID].FillTrigger(aTile->GetLocalTriggerPrimitive());;
3029           RootOutput->cd();
3030         }
3031 
3032         ithSpectra=hSpectra.find(currCellID);
3033         if (ithSpectra!=hSpectra.end()){
3034           ithSpectra->second.FillHGCROC(adc, toa, tot);
3035         } else {
3036           RootOutputHist->cd("IndividualCells");
3037           hSpectra[currCellID]=TileSpectra("mip1st",currCellID,calib.GetTileCalib(currCellID),event.GetROtype(),debug);
3038           hSpectra[currCellID].FillHGCROC(adc, toa, tot);
3039           RootOutput->cd();
3040         }
3041         
3042         // only fill tile spectra if 4 surrounding tiles on average are compatible with muon response
3043         if (localMuonTrigg){
3044           aTile->SetLocalTriggerBit(1);
3045           ithSpectraTrigg=hSpectraTrigg.find(currCellID);
3046           ithSpectraTrigg->second.FillHGCROC(adc, toa, tot);
3047         }
3048       }
3049     }
3050     RootOutput->cd();
3051     TdataOut->Fill();
3052   }
3053   TdataOut->Write();
3054   TsetupIn->CloneTree()->Write();
3055   
3056   //***********************************************************************************************
3057   //***** Monitoring histos for fits results of 2nd iteration ******************
3058   //***********************************************************************************************
3059   RootOutputHist->cd();
3060   int maxChannelPerLayer        = (setup->GetNMaxColumn()+1)*(setup->GetNMaxRow()+1);
3061   // monitoring trigger 
3062   TH2D* hmipTriggers              = new TH2D( "hmipTriggers","muon triggers; layer; brd channel; counts ",
3063                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
3064   hmipTriggers->SetDirectory(0);
3065   TH2D* hSuppresionNoise          = new TH2D( "hSuppresionNoise","S/B noise region; layer; brd channel; S/B noise region",
3066                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
3067   hSuppresionNoise->SetDirectory(0);
3068   TH2D* hSuppresionSignal         = new TH2D( "hSuppresionSignal","S/B signal region; layer; brd channel; S/B signal region",
3069                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
3070   hSuppresionSignal->SetDirectory(0);
3071 
3072   // monitoring 2nd iteration mip fits
3073   TH2D* hspectraHGMaxVsLayer   = new TH2D( "hspectraHGMaxVsLayer","Max High Gain; layer; brd channel; Max_{HG, 2^{nd}} (arb. units) ",
3074                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
3075   hspectraHGMaxVsLayer->SetDirectory(0);
3076   TH2D* hspectraHGFWHMVsLayer   = new TH2D( "hspectraHGFWHMVsLayer","FWHM High Gain; layer; brd channel; FWHM_{HG, 2^{nd}} (arb. units) ",
3077                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
3078   hspectraHGFWHMVsLayer->SetDirectory(0);
3079   TH2D* hspectraHGLMPVVsLayer   = new TH2D( "hspectraHGLMPVVsLayer","MPV High Gain; layer; brd channel; MPV_{HG, 2^{nd}} (arb. units) ",
3080                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
3081   hspectraHGLMPVVsLayer->SetDirectory(0);
3082   TH2D* hspectraHGLSigmaVsLayer = new TH2D( "hspectraHGLSigmaVsLayer","Sigma Landau High Gain; layer; brd channel; #sigma_{L,HG, 2^{nd}} (arb. units) ",
3083                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
3084   hspectraHGLSigmaVsLayer->SetDirectory(0);
3085   TH2D* hspectraHGGSigmaVsLayer = new TH2D( "hspectraHGGSigmaVsLayer","Sigma Gauss High Gain; layer; brd channel; #sigma_{G,HG, 2^{nd}} (arb. units) ",
3086                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
3087   hspectraHGGSigmaVsLayer->SetDirectory(0);
3088   TH2D* hspectraLGMaxVsLayer   = new TH2D( "hspectraLGMaxVsLayer","Max High Gain; layer; brd channel; Max_{LG, 2^{nd}} (arb. units) ",
3089                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
3090   hspectraLGMaxVsLayer->SetDirectory(0);
3091   TH2D* hspectraLGFWHMVsLayer   = new TH2D( "hspectraLGFWHMVsLayer","FWHM High Gain; layer; brd channel; FWHM_{LG, 2^{nd}} (arb. units) ",
3092                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
3093   hspectraLGFWHMVsLayer->SetDirectory(0);
3094   TH2D* hspectraLGLMPVVsLayer   = new TH2D( "hspectraLGLMPVVsLayer","MPV High Gain; layer; brd channel; MPV_{LG, 2^{nd}} (arb. units) ",
3095                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
3096   hspectraLGLMPVVsLayer->SetDirectory(0);
3097   TH2D* hspectraLGLSigmaVsLayer = new TH2D( "hspectraLGLSigmaVsLayer","Sigma Landau High Gain; layer; brd channel; #sigma_{L,LG, 2^{nd}} (arb. units) ",
3098                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
3099   hspectraLGLSigmaVsLayer->SetDirectory(0);
3100   TH2D* hspectraLGGSigmaVsLayer = new TH2D( "hspectraLGGSigmaVsLayer","Sigma Gauss High Gain; layer; brd channel; #sigma_{G,LG, 2^{nd}} (arb. units) ",
3101                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
3102   hspectraLGGSigmaVsLayer->SetDirectory(0);
3103 
3104   TH1D* hMaxHG             = new TH1D( "hMaxHG","Max High Gain ;Max_{HG} (arb. units) ; counts ",
3105                                             2000, -0.5, 2000-0.5);
3106   hMaxHG->SetDirectory(0);
3107   TH1D* hMaxLG             = new TH1D( "hMaxLG","Max Low Gain ;Max_{LG} (arb. units) ; counts ",
3108                                             400, -0.5, 400-0.5);
3109   hMaxLG->SetDirectory(0);
3110 
3111 
3112   int currCells = 0;
3113   double* parameters    = new double[6];
3114   double* parErrAndRes  = new double[6];
3115   bool isGood;
3116   double meanSB_NoiseR  = 0;
3117   double meanSB_SigR    = 0;
3118   if ( debug > 0)
3119     std::cout << "============================== start fitting improved iteration" << std::endl;  
3120   
3121   for(ithSpectraTrigg=hSpectraTrigg.begin(); ithSpectraTrigg!=hSpectraTrigg.end(); ++ithSpectraTrigg){
3122     if (currCells%20 == 0 && currCells > 0 && debug > 0)
3123       std::cout << "============================== cell " <<  currCells << " / " << hSpectraTrigg.size() << " cells" << std::endl;
3124     currCells++;
3125     for (int p = 0; p < 6; p++){
3126       parameters[p]   = 0;
3127       parErrAndRes[p] = 0;
3128     }
3129     isGood=ithSpectraTrigg->second.FitMipHG(parameters, parErrAndRes, debug, yearData, true, calib.GetVov(), averageScale);
3130     
3131     long cellID     = ithSpectraTrigg->second.GetCellID();
3132     int layer       = setup->GetLayer(cellID);
3133     int chInLayer   = setup->GetChannelInLayer(cellID);    
3134     int bin2D       = hspectraHGMaxVsLayer->FindBin(layer,chInLayer);
3135 
3136     double pedSigHG = 0;
3137     double maxBin   = 0;
3138     if (typeRO == ReadOut::Type::Caen){
3139       pedSigHG = calib.GetPedestalSigH(cellID);
3140       maxBin   = 3800;
3141     } else {
3142       pedSigHG = 20;
3143       maxBin   = 1024;
3144     }
3145 
3146     Int_t binNLow   = ithSpectraTrigg->second.GetHG()->FindBin(-1*pedSigHG);
3147     Int_t binNHigh  = ithSpectraTrigg->second.GetHG()->FindBin(3*pedSigHG);
3148     Int_t binSHigh  = ithSpectraTrigg->second.GetHG()->FindBin(maxBin);
3149     
3150     double S_NoiseR = ithSpectraTrigg->second.GetHG()->Integral(binNLow, binNHigh);
3151     double S_SigR   = ithSpectraTrigg->second.GetHG()->Integral(binNHigh, binSHigh);
3152     
3153     ithSpectra      = hSpectra.find(cellID);
3154     double B_NoiseR = ithSpectra->second.GetHG()->Integral(binNLow , binNHigh);
3155     double B_SigR   = ithSpectra->second.GetHG()->Integral(binNHigh, binSHigh);
3156     
3157     double SB_NoiseR  = (B_NoiseR != 0.) ? S_NoiseR/B_NoiseR : 0;
3158     double SB_SigR    = (B_SigR != 0.) ? S_SigR/B_SigR : 0;
3159     
3160     meanSB_NoiseR += SB_NoiseR;
3161     meanSB_SigR += SB_SigR;
3162     
3163     hmipTriggers->SetBinContent(bin2D, ithSpectraTrigg->second.GetHG()->GetEntries());
3164     hSuppresionNoise->SetBinContent(bin2D, SB_NoiseR);
3165     hSuppresionSignal->SetBinContent(bin2D, SB_SigR);
3166     if (isGood){
3167       hspectraHGMaxVsLayer->SetBinContent(bin2D, parameters[4]);
3168       hspectraHGFWHMVsLayer->SetBinContent(bin2D, parameters[5]);
3169       hspectraHGLMPVVsLayer->SetBinContent(bin2D, parameters[1]);
3170       hspectraHGLMPVVsLayer->SetBinError(bin2D, parErrAndRes[1]);
3171       hspectraHGLSigmaVsLayer->SetBinContent(bin2D, parameters[0]);
3172       hspectraHGLSigmaVsLayer->SetBinError(bin2D, parErrAndRes[0]);
3173       hspectraHGGSigmaVsLayer->SetBinContent(bin2D, parameters[3]);
3174       hspectraHGGSigmaVsLayer->SetBinError(bin2D, parErrAndRes[3]);
3175       hMaxHG->Fill(parameters[4]);
3176     }
3177     
3178     if (typeRO == ReadOut::Type::Caen) {
3179       for (int p = 0; p < 6; p++){
3180         parameters[p]   = 0;
3181         parErrAndRes[p] = 0;
3182       }
3183       isGood=ithSpectraTrigg->second.FitMipLG(parameters, parErrAndRes, debug, yearData, true, averageScaleLow);
3184       if (isGood){
3185         hspectraLGMaxVsLayer->SetBinContent(bin2D, parameters[4]);
3186         hspectraLGFWHMVsLayer->SetBinContent(bin2D, parameters[5]);
3187         hspectraLGLMPVVsLayer->SetBinContent(bin2D, parameters[1]);
3188         hspectraLGLMPVVsLayer->SetBinError(bin2D, parErrAndRes[1]);
3189         hspectraLGLSigmaVsLayer->SetBinContent(bin2D, parameters[0]);
3190         hspectraLGLSigmaVsLayer->SetBinError(bin2D, parErrAndRes[0]);
3191         hspectraLGGSigmaVsLayer->SetBinContent(bin2D, parameters[3]);
3192         hspectraLGGSigmaVsLayer->SetBinError(bin2D, parErrAndRes[3]);
3193         hMaxLG->Fill(parameters[4]);
3194       }
3195     }
3196   }
3197   if ( debug > 0)
3198     std::cout << "============================== done fitting improved iteration" << std::endl;
3199 
3200   
3201   meanSB_NoiseR = meanSB_NoiseR/hSpectraTrigg.size();
3202   meanSB_SigR   = meanSB_SigR/hSpectraTrigg.size();
3203   
3204   RootOutput->cd();
3205   
3206   if (IsCalibSaveToFile()){
3207     TString fileCalibPrint = RootOutputName;
3208     fileCalibPrint         = fileCalibPrint.ReplaceAll(".root","_calib.txt");
3209     calib.PrintCalibToFile(fileCalibPrint);
3210   }
3211   
3212   TcalibOut->Fill();
3213   TcalibOut->Write();
3214   int actChA                     = 0;
3215   double averageScaleUpdated     = calib.GetAverageScaleHigh(actChA);
3216   double averageScaleUpdatedLow  = 0.;
3217   double meanFWHMHG              = calib.GetAverageScaleWidthHigh();
3218   double meanFWHMLG              = 0.;
3219 
3220   if (typeRO == ReadOut::Type::Caen) {
3221     averageScaleUpdatedLow  = calib.GetAverageScaleLow();
3222     meanFWHMLG              = calib.GetAverageScaleWidthLow();
3223   }
3224   
3225   std::cout << "average input HG mip: " << averageScale << " LG mip: "<< averageScaleLow << "\t act. ch: "<< actChI<< std::endl;
3226   std::cout << "average updated HG mip: " << averageScaleUpdated << " LG mip: "<< averageScaleUpdatedLow << "\t act. ch: "<< actChA<< std::endl;
3227 
3228   RootOutput->Close();
3229 
3230 
3231   RootOutputHist->cd("IndividualCellsTrigg");
3232     for(ithSpectra=hSpectraTrigg.begin(); ithSpectra!=hSpectraTrigg.end(); ++ithSpectra){
3233       ithSpectra->second.Write(true);
3234     }
3235   RootOutputHist->cd();
3236     
3237     hspectraHGMaxVsLayer->Write();
3238     hspectraHGFWHMVsLayer->Write();
3239     hspectraHGLMPVVsLayer->Write();
3240     hspectraHGLSigmaVsLayer->Write();
3241     hspectraHGGSigmaVsLayer->Write();
3242     hMaxHG->Write();
3243     
3244     if (typeRO == ReadOut::Type::Caen){
3245       hspectraLGMaxVsLayer->Write();
3246       hspectraLGFWHMVsLayer->Write();
3247       hspectraLGLMPVVsLayer->Write();
3248       hspectraLGLSigmaVsLayer->Write();
3249       hspectraLGGSigmaVsLayer->Write();
3250       hMaxLG->Write();
3251     }
3252     hmipTriggers->Write();
3253     hSuppresionNoise->Write();
3254     hSuppresionSignal->Write();
3255   // fill calib tree & write it
3256   // close open root files
3257   RootOutputHist->Write();
3258   RootOutputHist->Close();
3259 
3260   RootInput->Close();
3261 
3262   // Get run info object
3263   std::map<int,RunInfo>::iterator it=ri.find(runNr);
3264   
3265   // create directory for plot output
3266   TString outputDirPlots = GetPlotOutputDir();
3267   gSystem->Exec("mkdir -p "+outputDirPlots);
3268   
3269   //**********************************************************************
3270   // Create canvases for channel overview plotting
3271   //**********************************************************************
3272   Double_t textSizeRel = 0.035;
3273   StyleSettingsBasics("pdf");
3274   SetPlotStyle();
3275 
3276   TCanvas* canvas2DCorr = new TCanvas("canvasCorrPlots","",0,0,1450,1300);  // gives the page size
3277   DefaultCancasSettings( canvas2DCorr, 0.08, 0.13, 0.045, 0.07);
3278 
3279   canvas2DCorr->SetLogz(0);
3280   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) );
3281   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));
3282   PlotSimple2D( canvas2DCorr, hspectraHGLMPVVsLayer, -10000, -10000, textSizeRel, Form("%s/HG_LandMPVMip.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
3283   PlotSimple2D( canvas2DCorr, hspectraHGLSigmaVsLayer, -10000, -10000, textSizeRel, Form("%s/HG_LandSigMip.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
3284   PlotSimple2D( canvas2DCorr, hspectraHGGSigmaVsLayer, -10000, -10000, textSizeRel, Form("%s/HG_GaussSigMip.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
3285   canvas2DCorr->SetLogz(1);
3286   TString drawOpt = "colz";
3287   if (yearData == 2023)
3288     drawOpt = "colz,text";
3289   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));
3290   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));
3291   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));
3292   
3293   if (typeRO == ReadOut::Type::Caen){
3294     canvas2DCorr->SetLogz(0);
3295     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));
3296     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));
3297     PlotSimple2D( canvas2DCorr, hspectraLGLMPVVsLayer, -10000, -10000, textSizeRel, Form("%s/LG_LandMPVMip.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
3298     PlotSimple2D( canvas2DCorr, hspectraLGLSigmaVsLayer, -10000, -10000, textSizeRel, Form("%s/LG_LandSigMip.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
3299     PlotSimple2D( canvas2DCorr, hspectraLGGSigmaVsLayer, -10000, -10000, textSizeRel, Form("%s/LG_GaussSigMip.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
3300   }
3301   //***********************************************************************************************************
3302   //********************************* 8 Panel overview plot  **************************************************
3303   //***********************************************************************************************************
3304   //*****************************************************************
3305     // Test beam geometry (beam coming from viewer)
3306     //===========================================================
3307     //||    8 (4)    ||    7 (5)   ||    6 (6)   ||    5 (7)   ||  row 0
3308     //===========================================================
3309     //||    1 (0)    ||    2 (1)   ||    3 (2)   ||    4 (3)   ||  row 1
3310     //===========================================================
3311     //    col 0     col 1       col 2     col  3
3312     // rebuild pad geom in similar way (numbering -1)
3313   //*****************************************************************
3314   TCanvas* canvas8Panel;
3315   TPad* pad8Panel[8];
3316   Double_t topRCornerX[8];
3317   Double_t topRCornerY[8];
3318   Int_t textSizePixel = 30;
3319   Double_t relSize8P[8];
3320   CreateCanvasAndPadsFor8PannelTBPlot(canvas8Panel, pad8Panel,  topRCornerX, topRCornerY, relSize8P, textSizePixel);
3321  
3322   calib.PrintGlobalInfo();  
3323   Double_t maxHG = ReturnMipPlotRangeDepVov(calib.GetVov(),true, typeRO);
3324   Double_t maxLG = ReturnMipPlotRangeDepVov(calib.GetVov(),false, typeRO);
3325   std::cout << "plotting single layers" << std::endl;
3326   for (Int_t l = 0; l < setup->GetNMaxLayer()+1; l++){   
3327     for (Int_t m = 0; m < setup->GetNMaxModule()+1; m++){ 
3328       if (l%10 == 0 && l > 0 && debug > 0)
3329         std::cout << "============================== layer " <<  l << " / " << setup->GetNMaxLayer() << " layers" << std::endl;
3330       PlotMipWithFits8MLayer (canvas8Panel,pad8Panel, topRCornerX, topRCornerY, relSize8P, textSizePixel, 
3331                                 hSpectra, hSpectraTrigg, setup, true, -100, maxHG, 1.2, l, m,
3332                                 Form("%s/MIP_HG_Mod%02d_Layer%02d.%s" ,outputDirPlots.Data(), m, l, plotSuffix.Data()), it->second);
3333       Double_t maxTriggPPlot = maxHG*2;
3334       if (typeRO != ReadOut::Type::Caen)
3335         maxTriggPPlot = 500;
3336       PlotTriggerPrimWithFits8MLayer (canvas8Panel,pad8Panel, topRCornerX, topRCornerY, relSize8P, textSizePixel, 
3337                                         hSpectraTrigg, setup, averageScale, factorMinTrigg, factorMaxTrigg,
3338                                         0, maxTriggPPlot, 1.2, l, m, Form("%s/TriggPrimitive_Mod%02d_Layer%02d.%s" ,outputDirPlots.Data(), m, l, plotSuffix.Data()), it->second);
3339       if (typeRO == ReadOut::Type::Caen){
3340         PlotMipWithFits8MLayer (canvas8Panel,pad8Panel, topRCornerX, topRCornerY, relSize8P, textSizePixel, 
3341                                   hSpectra, hSpectraTrigg, setup, false, -20, maxLG, 1.2, l, m,
3342                                   Form("%s/MIP_LG_Mod%02d_Layer%02d.%s" ,outputDirPlots.Data(), m, l, plotSuffix.Data()), it->second);
3343       }
3344     }
3345   }
3346   std::cout << "done plotting" << std::endl;
3347   
3348   if (ExtPlot > 0){
3349     TString outputDirPlotsSingle = Form("%s/SingleTiles",outputDirPlots.Data());
3350     gSystem->Exec("mkdir -p "+outputDirPlotsSingle);
3351 
3352     
3353     TCanvas* canvasSTile = new TCanvas("canvasSignleTile","",0,0,1600,1300);  // gives the page size
3354     DefaultCancasSettings( canvasSTile, 0.08, 0.01, 0.01, 0.082);
3355 
3356     int counter = 0;
3357     for(ithSpectraTrigg=hSpectraTrigg.begin(); ithSpectraTrigg!=hSpectraTrigg.end(); ++ithSpectraTrigg){
3358       counter++;
3359       long cellID     = ithSpectraTrigg->second.GetCellID();
3360       int row = setup->GetRow(cellID);
3361       int col = setup->GetColumn(cellID);
3362       int lay = setup->GetLayer(cellID);
3363       int mod = setup->GetModule(cellID);
3364 
3365       PlotMipWithFitsSingleTile (canvasSTile, 0.95,  0.95, Double_t(textSizePixel)*2/1600., textSizePixel*2, 
3366                                 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);
3367     }
3368   }
3369   return true;
3370 }
3371 
3372 
3373 
3374 //***********************************************************************************************
3375 //*********************** strip noise from sample and fit ***************************************
3376 //***********************************************************************************************
3377 bool Analyses::GetNoiseSampleAndRefitPedestal(void){
3378   std::cout<<"GetNoiseSampleAndRefitPedestal"<<std::endl;
3379   
3380   std::map<int,RunInfo> ri=readRunInfosFromFile(RunListInputName.Data(),debug,0);
3381   
3382   std::map<int,TileSpectra> hSpectra;
3383   std::map<int,TileSpectra> hSpectraTrigg;
3384   std::map<int, TileSpectra>::iterator ithSpectra;
3385   std::map<int, TileSpectra>::iterator ithSpectraTrigg;
3386   
3387   std::cout << "Additional Output with histos being created: " << RootOutputNameHist.Data() << std::endl;
3388   if(Overwrite){
3389     std::cout << "recreating file with hists" << std::endl;
3390     RootOutputHist = new TFile(RootOutputNameHist.Data(),"RECREATE");
3391   } else{
3392     std::cout << "newly creating file with hists" << std::endl;
3393     RootOutputHist = new TFile(RootOutputNameHist.Data(),"CREATE");
3394   }
3395     
3396   // setup trigger sel
3397   double factorMinTrigg   = 0.5;
3398   if(yearData == 2023)
3399     factorMinTrigg        = 0.1;
3400   // create HG and LG histo's per channel
3401   TH2D* hspectraHGvsCellID      = new TH2D( "hNoiseTriggeredSpectraHG_vsCellID","Noise trigg ADC spectrum High Gain vs CellID; cell ID; ADC_{HG} (arb. units) ",
3402                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5, 4000,0,4000);
3403   hspectraHGvsCellID->SetDirectory(0);
3404   TH2D* hspectraLGvsCellID      = new TH2D( "hNoiseTriggeredSpectraLG_vsCellID","Noise trigg ADC spectrum Low Gain vs CellID; cell ID; ADC_{LG} (arb. units)  ",
3405                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5, 4000,0,4000);
3406   hspectraLGvsCellID->SetDirectory(0);
3407 
3408   
3409   RootOutputHist->mkdir("IndividualCells");
3410   RootOutputHist->mkdir("IndividualCellsTrigg");
3411   RootOutputHist->cd("IndividualCellsTrigg");  
3412   
3413   //***********************************************************************************************
3414   //************************* first pass over tree to extract spectra *****************************
3415   //***********************************************************************************************  
3416   TcalibIn->GetEntry(0);
3417   // check whether calib should be overwritten based on external text file
3418   if (OverWriteCalib){
3419     calib.ReadCalibFromTextFile(ExternalCalibFile,debug);
3420   }
3421   
3422   int evts=TdataIn->GetEntries();
3423   int runNr = -1;
3424   int actCh = 0;
3425   double averageScale     = calib.GetAverageScaleHigh(actCh);
3426   double averageScaleLow  = calib.GetAverageScaleLow();
3427   std::cout << "average HG mip: " << averageScale << " LG mip: "<< averageScaleLow << std::endl;
3428   for(int i=0; i<evts; i++){
3429     TdataIn->GetEntry(i);
3430     if (i == 0)runNr = event.GetRunNumber();
3431     TdataIn->GetEntry(i);
3432     if (i%5000 == 0 && i > 0 && debug > 0) std::cout << "Reading " <<  i << " / " << evts << " events" << std::endl;
3433     for(int j=0; j<event.GetNTiles(); j++){
3434       Caen* aTile=(Caen*)event.GetTile(j);
3435       if (i == 0 && debug > 2) std::cout << ((TString)setup->DecodeCellID(aTile->GetCellID())).Data() << std::endl;
3436       long currCellID = aTile->GetCellID();
3437       
3438       // read current tile
3439       ithSpectraTrigg=hSpectraTrigg.find(aTile->GetCellID());
3440       // estimate local muon trigger
3441       bool localNoiseTrigg = event.InspectIfNoiseTrigg(currCellID, averageScale, factorMinTrigg);
3442       
3443       if(ithSpectraTrigg!=hSpectraTrigg.end()){
3444         ithSpectraTrigg->second.FillTrigger(aTile->GetLocalTriggerPrimitive());
3445       } else {
3446         RootOutputHist->cd("IndividualCellsTrigg");
3447         hSpectraTrigg[currCellID]=TileSpectra("mipTrigg",currCellID,calib.GetTileCalib(currCellID),event.GetROtype(),debug);
3448         hSpectraTrigg[currCellID].FillTrigger(aTile->GetLocalTriggerPrimitive());;
3449         RootOutput->cd();
3450       }
3451       
3452       ithSpectra=hSpectra.find(aTile->GetCellID());
3453       if (ithSpectra!=hSpectra.end()){
3454         ithSpectra->second.FillSpectraCAEN(aTile->GetADCLow(),aTile->GetADCHigh());
3455       } else {
3456         RootOutputHist->cd("IndividualCells");
3457         hSpectra[currCellID]=TileSpectra("mip1st",currCellID,calib.GetTileCalib(currCellID),event.GetROtype(),debug);
3458         hSpectra[aTile->GetCellID()].FillSpectraCAEN(aTile->GetADCLow(),aTile->GetADCHigh());;
3459 
3460         RootOutput->cd();
3461       }
3462      
3463       // only fill tile spectra if X surrounding tiles on average are compatible with pure noise
3464       if (localNoiseTrigg){
3465         aTile->SetLocalTriggerBit(2);
3466         ithSpectraTrigg=hSpectraTrigg.find(aTile->GetCellID());
3467         ithSpectraTrigg->second.FillSpectraCAEN(aTile->GetADCLow(),aTile->GetADCHigh());
3468         
3469         hspectraHGvsCellID->Fill(aTile->GetCellID(), aTile->GetADCHigh());
3470         hspectraLGvsCellID->Fill(aTile->GetCellID(), aTile->GetADCLow());
3471       }
3472     }
3473     RootOutput->cd();
3474     TdataOut->Fill();
3475   }
3476   TdataOut->Write();
3477   TsetupIn->CloneTree()->Write();
3478   
3479   //***********************************************************************************************
3480   //***** Monitoring histos for fits results of 2nd iteration ******************
3481   //***********************************************************************************************
3482   RootOutputHist->cd();
3483   int maxChannelPerLayer        = (setup->GetNMaxColumn()+1)*(setup->GetNMaxRow()+1);
3484   // monitoring trigger 
3485   TH2D* hnoiseTriggers              = new TH2D( "hnoiseTriggers","muon triggers; layer; brd channel; counts ",
3486                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
3487   hnoiseTriggers->SetDirectory(0);
3488   TH1D* hMeanPedHGvsCellID      = new TH1D( "hMeanPedHG_vsCellID","mean Ped High Gain vs CellID ; cell ID; #mu_{noise, HG} (arb. units) ",
3489                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5);
3490   hMeanPedHGvsCellID->SetDirectory(0);
3491   TH2D* hspectraHGMeanVsLayer   = new TH2D( "hspectraHGMeanVsLayer","Mean Ped High Gain vs CellID; layer; brd channel; #mu_{Ped HG} (arb. units) ",
3492                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
3493   hspectraHGMeanVsLayer->SetDirectory(0);
3494   TH2D* hspectraHGSigmaVsLayer  = new TH2D( "hspectraHGSigmaVsLayer","Mean Ped High Gain vs CellID; layer; brd channel; #sigma_{Ped HG} (arb. units) ",
3495                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
3496   hspectraHGSigmaVsLayer->SetDirectory(0);
3497   TH1D* hMeanPedLGvsCellID      = new TH1D( "hMeanPedLG_vsCellID","mean Ped Low Gain vs CellID ; cell ID; #mu_{noise, LG} (arb. units) ",
3498                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5);
3499   hMeanPedLGvsCellID->SetDirectory(0);
3500   TH2D* hspectraLGMeanVsLayer   = new TH2D( "hspectraLGMeanVsLayer","Mean Ped Low Gain vs CellID; layer; brd channel; #mu_{PED LG} (arb. units) ",
3501                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
3502   hspectraLGMeanVsLayer->SetDirectory(0);
3503   TH2D* hspectraLGSigmaVsLayer  = new TH2D( "hspectraLGSigmaVsLayer","Mean Ped Low Gain vs CellID; layer; brd channel; #sigma_{Ped LG} (arb. units)",
3504                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
3505   hspectraLGSigmaVsLayer->SetDirectory(0);
3506 
3507   if ( debug > 0)
3508     std::cout << "============================== starting fitting" << std::endl;
3509 
3510   int currCells = 0;
3511   double* parameters    = new double[6];
3512   for(ithSpectraTrigg=hSpectraTrigg.begin(); ithSpectraTrigg!=hSpectraTrigg.end(); ++ithSpectraTrigg){
3513     for (int p = 0; p < 6; p++){
3514       parameters[p]   = 0;
3515     }
3516 
3517     if (currCells%20 == 0 && currCells > 0 && debug > 0)
3518       std::cout << "============================== cell " <<  currCells << " / " << hSpectraTrigg.size() << " cells" << std::endl;
3519     currCells++;
3520     if ( debug > 2) std::cout << ((TString)setup->DecodeCellID(ithSpectraTrigg->second.GetCellID())).Data() << std::endl;
3521     ithSpectraTrigg->second.FitNoise(parameters, yearData, true);
3522     hMeanPedHGvsCellID->SetBinContent(hMeanPedHGvsCellID->FindBin(ithSpectraTrigg->second.GetCellID()), parameters[4]);
3523     hMeanPedHGvsCellID->SetBinError  (hMeanPedHGvsCellID->FindBin(ithSpectraTrigg->second.GetCellID()), parameters[6]);
3524     hMeanPedLGvsCellID->SetBinContent(hMeanPedLGvsCellID->FindBin(ithSpectraTrigg->second.GetCellID()), parameters[0]);
3525     hMeanPedLGvsCellID->SetBinError  (hMeanPedLGvsCellID->FindBin(ithSpectraTrigg->second.GetCellID()), parameters[2]);
3526     
3527     int layer     = setup->GetLayer(ithSpectraTrigg->second.GetCellID());
3528     int chInLayer = setup->GetChannelInLayer(ithSpectraTrigg->second.GetCellID());
3529   
3530     hspectraHGMeanVsLayer->SetBinContent(hspectraHGMeanVsLayer->FindBin(layer,chInLayer), parameters[4]);
3531     hspectraHGMeanVsLayer->SetBinError(hspectraHGMeanVsLayer->FindBin(layer,chInLayer), parameters[5]);
3532     hspectraHGSigmaVsLayer->SetBinContent(hspectraHGSigmaVsLayer->FindBin(layer,chInLayer), parameters[6]);
3533     hspectraHGSigmaVsLayer->SetBinError(hspectraHGSigmaVsLayer->FindBin(layer,chInLayer), parameters[7]);
3534     hspectraLGMeanVsLayer->SetBinContent(hspectraLGMeanVsLayer->FindBin(layer,chInLayer), parameters[0]);
3535     hspectraLGMeanVsLayer->SetBinError(hspectraLGMeanVsLayer->FindBin(layer,chInLayer), parameters[1]);
3536     hspectraLGSigmaVsLayer->SetBinContent(hspectraLGSigmaVsLayer->FindBin(layer,chInLayer), parameters[2]);
3537     hspectraLGSigmaVsLayer->SetBinError(hspectraLGSigmaVsLayer->FindBin(layer,chInLayer), parameters[3]);
3538     
3539     hnoiseTriggers->SetBinContent(hnoiseTriggers->FindBin(layer,chInLayer), ithSpectraTrigg->second.GetHG()->GetEntries());
3540   }
3541   if ( debug > 0)
3542     std::cout << "============================== done fitting" << std::endl;
3543   
3544   RootOutput->cd();
3545   
3546   if (IsCalibSaveToFile()){
3547     TString fileCalibPrint = RootOutputName;
3548     fileCalibPrint         = fileCalibPrint.ReplaceAll(".root","_calib.txt");
3549     calib.PrintCalibToFile(fileCalibPrint);
3550   }
3551 
3552 
3553   TcalibOut->Fill();
3554   TcalibOut->Write();
3555   
3556   RootOutput->Write();
3557   RootOutput->Close();
3558 
3559 
3560   RootOutputHist->cd("IndividualCellsTrigg");
3561     for(ithSpectra=hSpectraTrigg.begin(); ithSpectra!=hSpectraTrigg.end(); ++ithSpectra){
3562       ithSpectra->second.Write(true);
3563     }
3564   RootOutputHist->cd();
3565     
3566     hMeanPedHGvsCellID->Write();
3567     hMeanPedLGvsCellID->Write();
3568     hspectraHGMeanVsLayer->Write();
3569     hspectraHGSigmaVsLayer->Write();
3570     hspectraLGMeanVsLayer->Write(); 
3571     hspectraLGSigmaVsLayer->Write();
3572     hspectraHGvsCellID->Write();
3573     hspectraLGvsCellID->Write();
3574         
3575     hnoiseTriggers->Write();
3576   // close open root files
3577   RootOutputHist->Write();
3578   RootOutputHist->Close();
3579 
3580   RootInput->Close();
3581 
3582   // Get run info object
3583   std::map<int,RunInfo>::iterator it=ri.find(runNr);
3584   
3585   // create directory for plot output
3586   TString outputDirPlots = GetPlotOutputDir();
3587   gSystem->Exec("mkdir -p "+outputDirPlots);
3588   
3589   //**********************************************************************
3590   // Create canvases for channel overview plotting
3591   //**********************************************************************
3592   Double_t textSizeRel = 0.035;
3593   StyleSettingsBasics("pdf");
3594   SetPlotStyle();
3595 
3596   TCanvas* canvas2DCorr = new TCanvas("canvasCorrPlots","",0,0,1450,1300);  // gives the page size
3597   DefaultCancasSettings( canvas2DCorr, 0.08, 0.13, 0.045, 0.07);
3598 
3599   canvas2DCorr->SetLogz(0);
3600   PlotSimple2D( canvas2DCorr, hspectraHGMeanVsLayer, -10000, -10000, textSizeRel, Form("%s/HG_NoiseMean.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true );
3601   PlotSimple2D( canvas2DCorr, hspectraHGSigmaVsLayer, -10000, -10000, textSizeRel, Form("%s/HG_NoiseSigma.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
3602   PlotSimple2D( canvas2DCorr, hspectraLGMeanVsLayer, -10000, -10000, textSizeRel, Form("%s/LG_NoiseMean.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
3603   PlotSimple2D( canvas2DCorr, hspectraLGSigmaVsLayer, -10000, -10000, textSizeRel, Form("%s/LG_NoiseSigma.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
3604   canvas2DCorr->SetLogz(1);
3605   PlotSimple2D( canvas2DCorr, hspectraHGvsCellID, -10000, -10000, textSizeRel, Form("%s/HG_Noise.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
3606   PlotSimple2D( canvas2DCorr, hspectraLGvsCellID, -10000, -10000, textSizeRel, Form("%s/LG_Noise.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
3607   
3608   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));
3609   //***********************************************************************************************************
3610   //********************************* 8 Panel overview plot  **************************************************
3611   //***********************************************************************************************************
3612   //*****************************************************************
3613     // Test beam geometry (beam coming from viewer)
3614     //===========================================================
3615     //||    8 (4)    ||    7 (5)   ||    6 (6)   ||    5 (7)   ||  row 0
3616     //===========================================================
3617     //||    1 (0)    ||    2 (1)   ||    3 (2)   ||    4 (3)   ||  row 1
3618     //===========================================================
3619     //    col 0     col 1       col 2     col  3
3620     // rebuild pad geom in similar way (numbering -1)
3621   //*****************************************************************
3622   TCanvas* canvas8Panel;
3623   TPad* pad8Panel[8];
3624   Double_t topRCornerX[8];
3625   Double_t topRCornerY[8];
3626   Int_t textSizePixel = 30;
3627   Double_t relSize8P[8];
3628   CreateCanvasAndPadsFor8PannelTBPlot(canvas8Panel, pad8Panel,  topRCornerX, topRCornerY, relSize8P, textSizePixel);
3629  
3630   for (Int_t l = 0; l < setup->GetNMaxLayer()+1; l++){
3631     for (Int_t m = 0; m < setup->GetNMaxModule()+1; m++){    
3632       PlotNoiseAdvWithFits8MLayer (canvas8Panel,pad8Panel, topRCornerX, topRCornerY, relSize8P, textSizePixel, 
3633                                       hSpectra, hSpectraTrigg, true, 0, 450, 1.2, l, m,
3634                                       Form("%s/NoiseTrigg_HG_Mod%02d_Layer%02d.%s" ,outputDirPlots.Data(), m, l, plotSuffix.Data()), it->second);
3635     }
3636   }
3637 
3638   
3639   return true;
3640 }
3641 
3642 
3643 //***********************************************************************************************
3644 //*********************** Evaluate local triggers only and store ********************************
3645 //***********************************************************************************************
3646 bool Analyses::RunEvalLocalTriggers(void){
3647   std::cout<<"EvalLocalTriggers"<<std::endl;
3648 
3649   RootOutput->cd();
3650   std::cout << "starting to run trigger eval: " << TcalibIn <<  "\t" << TcalibIn->GetEntry(0) << std::endl;
3651   TcalibIn->GetEntry(0);
3652   // check whether calib should be overwritten based on external text file
3653   if (OverWriteCalib){
3654     calib.ReadCalibFromTextFile(ExternalCalibFile,debug);
3655   }
3656   
3657   int actCh1st = 0;
3658   double averageScale = calib.GetAverageScaleHigh(actCh1st);
3659   double avLGHGCorr   = calib.GetAverageLGHGCorr();
3660   std::cout << "average HG mip: " << averageScale << "\t active ch: "<< actCh1st<< std::endl;
3661   
3662   // setup local trigger sel
3663   TRandom3* rand    = new TRandom3();
3664   Int_t localTriggerTiles = 4;
3665   if (yearData == 2023){
3666     localTriggerTiles = 6;
3667   }
3668   double factorMinTrigg   = 0.8;
3669   double factorMinTriggNoise = 0.2;
3670   double factorMaxTrigg   = 2.;
3671   if (yearData == 2023){
3672     factorMinTrigg    = 0.9;
3673     factorMaxTrigg    = 2.;
3674   }
3675   
3676   int outCount      = 1000;
3677   int evts=TdataIn->GetEntries();
3678   
3679   if (maxEvents == -1){
3680     maxEvents = evts;
3681   } else {
3682     std::cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl;
3683     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;
3684     std::cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl;
3685   }
3686   
3687   if (evts < 10000)
3688     outCount  = 500;
3689   if (evts > 100000)
3690     outCount  = 5000;
3691   int runNr = -1;
3692   for(int i=0; i<evts && i < maxEvents; i++){
3693     TdataIn->GetEntry(i);
3694     if(debug==1000){
3695         std::cerr<<event.GetTimeStamp()<<std::endl;
3696       }
3697     if (i == 0){
3698       runNr = event.GetRunNumber();
3699       std::cout<< "original run numbers calib: "<<calib.GetRunNumber() << "\t" << calib.GetRunNumberPed() << "\t" << calib.GetRunNumberMip() << std::endl;
3700       calib.SetRunNumber(runNr);
3701       calib.SetBeginRunTime(event.GetBeginRunTimeAlt());
3702       std::cout<< "reset run numbers calib: "<< calib.GetRunNumber() << "\t" << calib.GetRunNumberPed() << "\t" << calib.GetRunNumberMip() << std::endl;
3703     }
3704     
3705     if (i%outCount == 0 && debug > 0) std::cout << "Reading " <<  i << " / " << evts << " events" << std::endl;
3706     for(int j=0; j<event.GetNTiles(); j++){
3707       Caen* aTile=(Caen*)event.GetTile(j);      
3708       // calculate trigger primitives
3709       aTile->SetLocalTriggerPrimitive(event.CalculateLocalMuonTrigg(calib, rand, aTile->GetCellID(), localTriggerTiles, avLGHGCorr));
3710       bool localMuonTrigg   = event.InspectIfLocalMuonTrigg(aTile->GetCellID(), averageScale, factorMinTrigg, factorMaxTrigg);
3711       bool localNoiseTrigg  = event.InspectIfNoiseTrigg(aTile->GetCellID(), averageScale, factorMinTriggNoise);
3712       aTile->SetLocalTriggerBit(0);
3713       if (localMuonTrigg) aTile->SetLocalTriggerBit(1);
3714       if (localNoiseTrigg) aTile->SetLocalTriggerBit(2);
3715     }
3716     TdataOut->Fill();
3717   }
3718   TdataOut->Write();
3719   TsetupIn->CloneTree()->Write();
3720   
3721   if (IsCalibSaveToFile()){
3722     TString fileCalibPrint = RootOutputName;
3723     fileCalibPrint         = fileCalibPrint.ReplaceAll(".root","_calib.txt");
3724     calib.PrintCalibToFile(fileCalibPrint);
3725   }
3726   
3727   TcalibOut->Fill();
3728   TcalibOut->Write();
3729   
3730   RootOutput->Close();
3731   RootInput->Close();      
3732   
3733   return true;
3734 }
3735 
3736 //***********************************************************************************************
3737 //*********************** Calibration routine ***************************************************
3738 //***********************************************************************************************
3739 bool Analyses::Calibrate(void){
3740   std::cout<<"Calibrate"<<std::endl;
3741 
3742   std::map<int,RunInfo> ri=readRunInfosFromFile(RunListInputName.Data(),debug,0);
3743 
3744   std::cout << "Additional Output with histos being created: " << RootOutputNameHist.Data() << std::endl;
3745   if(Overwrite){
3746     std::cout << "recreating file with hists" << std::endl;
3747     RootOutputHist = new TFile(RootOutputNameHist.Data(),"RECREATE");
3748   } else{
3749     std::cout << "newly creating file with hists" << std::endl;
3750     RootOutputHist = new TFile(RootOutputNameHist.Data(),"CREATE");
3751   }
3752   
3753   // create HG and LG histo's per channel
3754   TH2D* hspectraHGCorrvsCellID      = new TH2D( "hspectraHGCorr_vsCellID","ADC spectrum High Gain corrected vs CellID; cell ID; ADC_{HG} (arb. units)  ; counts ",
3755                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5, 4000,-200,3800);
3756   hspectraHGCorrvsCellID->SetDirectory(0);
3757   TH2D* hspectraHGCorrvsCellIDNoise      = new TH2D( "hspectraHGCorr_vsCellID_Noise","ADC spectrum High Gain corrected vs CellID Noise; cell ID; ADC_{HG} (arb. units)  ; counts ",
3758                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5, 4000,-200,3800);
3759   hspectraHGCorrvsCellIDNoise->SetDirectory(0);
3760   TH2D* hspectraLGCorrvsCellID      = new TH2D( "hspectraLGCorr_vsCellID","ADC spectrum Low Gain corrected vs CellID; cell ID; ADC_{LG} (arb. units)  ; counts  ",
3761                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5, 4000,-200,3800);
3762   hspectraLGCorrvsCellID->SetDirectory(0);
3763   TH2D* hspectraLGCorrvsCellIDNoise      = new TH2D( "hspectraLGCorr_vsCellID_Noise","ADC spectrum Low Gain corrected vs CellID Noise; cell ID; ADC_{LG} (arb. units)  ; counts  ",
3764                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5, 4000,-200,3800);
3765   hspectraLGCorrvsCellIDNoise->SetDirectory(0);
3766   TH2D* hspectraHGvsCellID      = new TH2D( "hspectraHG_vsCellID","ADC spectrum High Gain vs CellID; cell ID; ADC_{HG} (arb. units)   ; counts ",
3767                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5, 4000,0,4000);
3768   hspectraHGvsCellID->SetDirectory(0);
3769   TH2D* hspectraLGvsCellID      = new TH2D( "hspectraLG_vsCellID","ADC spectrum Low Gain vs CellID; cell ID; ADC_{LG} (arb. units)  ; counts",
3770                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5, 4000,0,4000);
3771   hspectraLGvsCellID->SetDirectory(0);
3772   TH2D* hspectraEnergyvsCellID  = new TH2D( "hspectraEnergy_vsCellID","Energy vs CellID; cell ID; E (mip eq./tile)    ; counts",
3773                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5, 6000,0,200);
3774   hspectraEnergyvsCellID->SetDirectory(0);
3775   TH2D* hspectraEnergyTotvsNCells  = new TH2D( "hspectraTotEnergy_vsNCells","Energy vs CellID; N_{Cells}; E_{tot} (mip eq./tile) ; counts",
3776                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5, 6000,0,1000);
3777   hspectraEnergyTotvsNCells->SetDirectory(0);
3778 
3779   std::map<int,TileSpectra> hSpectra;
3780   std::map<int, TileSpectra>::iterator ithSpectra;
3781   std::map<int,TileSpectra> hSpectraNoise;
3782   std::map<int, TileSpectra>::iterator ithSpectraNoise;
3783   
3784   // entering histoOutput file
3785   RootOutputHist->mkdir("IndividualCells");
3786   RootOutputHist->cd("IndividualCells");
3787   RootOutputHist->mkdir("IndividualCellsNoise");
3788   RootOutputHist->cd("IndividualCellsNoise");
3789 
3790   Int_t runNr = -1;
3791   RootOutput->cd();
3792   std::cout << "starting to run calibration: " << TcalibIn <<  "\t" << TcalibIn->GetEntry(0) << std::endl;
3793   TcalibIn->GetEntry(0);
3794   // check whether calib should be overwritten based on external text file
3795   if (OverWriteCalib){
3796     calib.ReadCalibFromTextFile(ExternalCalibFile,debug);
3797   }
3798 
3799   int actCh1st = 0;
3800   double averageScale = calib.GetAverageScaleHigh(actCh1st);
3801   double avLGHGCorr   = calib.GetAverageLGHGCorr();
3802   std::cout << "average HG mip: " << averageScale << "\t active ch: "<< actCh1st<< std::endl;
3803   
3804   // setup local trigger sel
3805   TRandom3* rand    = new TRandom3();
3806   Int_t localTriggerTiles = 4;
3807   if (yearData == 2023){
3808     localTriggerTiles = 6;
3809   }
3810   double factorMinTrigg   = 0.8;
3811   double factorMinTriggNoise = 0.2;
3812   double factorMaxTrigg   = 2.;
3813   if (yearData == 2023){
3814     factorMinTrigg    = 0.9;
3815     factorMaxTrigg    = 2.;
3816   }
3817 
3818   
3819   double minMipFrac = 0.3;
3820   int corrHGADCSwap = 3500;
3821   int outCount      = 5000;
3822   int evts=TdataIn->GetEntries();
3823   
3824   if (maxEvents == -1){
3825     maxEvents = evts;
3826   } else {
3827     std::cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl;
3828     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;
3829     std::cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl;
3830   }
3831 
3832   if (evts < 10000)
3833     outCount  = 500;
3834   for(int i=0; i<evts && i < maxEvents; i++){
3835     if(debug==1000){
3836         std::cerr<<event.GetTimeStamp()<<std::endl;
3837       }
3838     TdataIn->GetEntry(i);
3839     if (i%outCount == 0 && debug > 0) std::cout << "Reading " <<  i << " / " << evts << " events" << std::endl;
3840     if (i == 0){
3841       runNr = event.GetRunNumber();
3842       std::cout<< "original run numbers calib: "<<calib.GetRunNumber() << "\t" << calib.GetRunNumberPed() << "\t" << calib.GetRunNumberMip() << std::endl;
3843       calib.SetRunNumber(runNr);
3844       calib.SetBeginRunTime(event.GetBeginRunTimeAlt());
3845       std::cout<< "reset run numbers calib: "<< calib.GetRunNumber() << "\t" << calib.GetRunNumberPed() << "\t" << calib.GetRunNumberMip() << std::endl;
3846     }
3847     double Etot = 0;
3848     int nCells  = 0;
3849     for(int j=0; j<event.GetNTiles(); j++){
3850       Caen* aTile=(Caen*)event.GetTile(j);
3851       // remove bad channels from output
3852       if (calib.GetBCCalib() && calib.GetBadChannel(aTile->GetCellID())!= 3 ){
3853         event.RemoveTile(aTile);
3854         j--;        
3855         continue;
3856       }
3857       
3858       double energy = 0;
3859       double corrHG = aTile->GetADCHigh()-calib.GetPedestalMeanH(aTile->GetCellID());
3860       double corrLG = aTile->GetADCLow()-calib.GetPedestalMeanL(aTile->GetCellID());
3861       double corrLG_HGeq = corrLG*calib.GetLGHGCorr(aTile->GetCellID()) + calib.GetLGHGCorrOff(aTile->GetCellID());
3862       if(corrHG<corrHGADCSwap){
3863         if(corrHG/calib.GetScaleHigh(aTile->GetCellID()) > minMipFrac){
3864           energy=corrHG/calib.GetScaleHigh(aTile->GetCellID());
3865         }
3866       } else {
3867         energy=corrLG/calib.GetCalcScaleLow(aTile->GetCellID());
3868       }
3869       if (debug > 1 && corrHG >= corrHGADCSwap-100 && corrHG < 4000-calib.GetPedestalMeanH(aTile->GetCellID())){
3870           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;
3871       }
3872       // calculate trigger primitives
3873       bool localMuonTrigg   = false;
3874       bool localNoiseTrigg  = false;
3875 
3876       if (!UseLocTriggFromFile()){
3877         aTile->SetLocalTriggerPrimitive(event.CalculateLocalMuonTrigg(calib, rand, aTile->GetCellID(), localTriggerTiles, avLGHGCorr));
3878         aTile->SetLocalTriggerBit(0);
3879         localMuonTrigg   = event.InspectIfLocalMuonTrigg(aTile->GetCellID(), averageScale, factorMinTrigg, factorMaxTrigg);
3880         localNoiseTrigg  = event.InspectIfNoiseTrigg(aTile->GetCellID(), averageScale, factorMinTriggNoise);
3881         if (localMuonTrigg) aTile->SetLocalTriggerBit(1);
3882         if (localNoiseTrigg) aTile->SetLocalTriggerBit(2);
3883       } else {
3884         if (aTile->GetLocalTriggerBit() == 2)  localNoiseTrigg = true;
3885       }
3886       
3887       hspectraHGvsCellID->Fill(aTile->GetCellID(), aTile->GetADCHigh());
3888       hspectraLGvsCellID->Fill(aTile->GetCellID(), aTile->GetADCLow());
3889       hspectraHGCorrvsCellID->Fill(aTile->GetCellID(), corrHG);
3890       hspectraLGCorrvsCellID->Fill(aTile->GetCellID(), corrLG);
3891       
3892       if (localNoiseTrigg){
3893         hspectraHGCorrvsCellIDNoise->Fill(aTile->GetCellID(), corrHG);
3894         hspectraLGCorrvsCellIDNoise->Fill(aTile->GetCellID(), corrLG);
3895       
3896       }
3897       ithSpectra=hSpectra.find(aTile->GetCellID());
3898       if(ithSpectra!=hSpectra.end()){
3899         ithSpectra->second.FillExtCAEN(corrLG,corrHG,energy,corrLG_HGeq);
3900       } else {
3901         RootOutputHist->cd("IndividualCells");
3902         hSpectra[aTile->GetCellID()]=TileSpectra("Calibrate",1,aTile->GetCellID(),calib.GetTileCalib(aTile->GetCellID()),event.GetROtype(),debug);
3903         hSpectra[aTile->GetCellID()].FillExtCAEN(corrLG,corrHG,energy,corrLG_HGeq);
3904         RootOutput->cd();
3905       }
3906 
3907       if (localNoiseTrigg){
3908         ithSpectraNoise=hSpectraNoise.find(aTile->GetCellID());
3909         if(ithSpectraNoise!=hSpectraNoise.end()){
3910           ithSpectraNoise->second.FillExtCAEN(corrLG,corrHG,energy,corrLG_HGeq);
3911         } else {
3912           RootOutputHist->cd("IndividualCellsNoise");
3913           hSpectraNoise[aTile->GetCellID()]=TileSpectra("CalibrateNoise",1,aTile->GetCellID(),calib.GetTileCalib(aTile->GetCellID()),event.GetROtype(),debug);
3914           hSpectraNoise[aTile->GetCellID()].FillExtCAEN(corrLG,corrHG,energy,corrLG_HGeq);
3915           RootOutput->cd();
3916         }
3917       }
3918       
3919       if(energy!=0){ 
3920         aTile->SetE(energy);
3921         hspectraEnergyvsCellID->Fill(aTile->GetCellID(), energy);
3922         Etot=Etot+energy;
3923         nCells++;
3924       } else {
3925         event.RemoveTile(aTile);
3926         j--;
3927       }
3928     }
3929     hspectraEnergyTotvsNCells->Fill(nCells,Etot);
3930     RootOutput->cd();
3931     TdataOut->Fill();
3932   }
3933   TdataOut->Write();
3934   TsetupIn->CloneTree()->Write();
3935   
3936   if (IsCalibSaveToFile()){
3937     TString fileCalibPrint = RootOutputName;
3938     fileCalibPrint         = fileCalibPrint.ReplaceAll(".root","_calib.txt");
3939     calib.PrintCalibToFile(fileCalibPrint);
3940   }
3941   
3942   TcalibOut->Fill();
3943   TcalibOut->Write();
3944   
3945   
3946   RootOutput->Close();
3947   RootInput->Close();      
3948   
3949   RootOutputHist->cd();
3950 
3951   hspectraHGvsCellID->Write();
3952   hspectraLGvsCellID->Write();
3953   hspectraHGCorrvsCellID->Write();
3954   hspectraLGCorrvsCellID->Write();
3955   hspectraHGCorrvsCellIDNoise->Write();
3956   hspectraLGCorrvsCellIDNoise->Write();
3957   hspectraEnergyvsCellID->Write();
3958   hspectraEnergyTotvsNCells->Write();
3959   
3960   TH1D* hspectraEnergyTot = (TH1D*)hspectraEnergyTotvsNCells->ProjectionY();
3961   hspectraEnergyTot->SetDirectory(0);
3962   TH1D* hspectraNCells = (TH1D*)hspectraEnergyTotvsNCells->ProjectionX();
3963   hspectraNCells->SetDirectory(0);
3964   hspectraEnergyTot->Write("hTotEnergy");
3965   hspectraNCells->Write("hNCells");
3966 
3967   RootOutputHist->cd("IndividualCells");
3968   for(ithSpectra=hSpectra.begin(); ithSpectra!=hSpectra.end(); ++ithSpectra){
3969     ithSpectra->second.FitLGHGCorr(debug,false);
3970     ithSpectra->second.WriteExt(true);
3971   }
3972   RootOutputHist->cd("IndividualCellsNoise");
3973   for(ithSpectraNoise=hSpectraNoise.begin(); ithSpectraNoise!=hSpectraNoise.end(); ++ithSpectraNoise){
3974     ithSpectraNoise->second.WriteExt(true);
3975   }
3976   
3977   RootOutputHist->Close();
3978   //**********************************************************************
3979   //********************* Plotting ***************************************
3980   //**********************************************************************
3981   // Get run info object
3982   std::map<int,RunInfo>::iterator it=ri.find(runNr);
3983   
3984   TString outputDirPlots = GetPlotOutputDir();
3985   gSystem->Exec("mkdir -p "+outputDirPlots);
3986   
3987   //**********************************************************************
3988   // Create canvases for channel overview plotting
3989   //**********************************************************************
3990   Double_t textSizeRel = 0.035;
3991   StyleSettingsBasics("pdf");
3992   SetPlotStyle();
3993   
3994   TCanvas* canvas2DCorr = new TCanvas("canvasCorrPlots","",0,0,1450,1300);  // gives the page size
3995   DefaultCancasSettings( canvas2DCorr, 0.08, 0.13, 0.045, 0.07);
3996   canvas2DCorr->SetLogz(1);
3997   PlotSimple2D( canvas2DCorr, hspectraHGvsCellID, -10000, -10000, textSizeRel, Form("%s/HG.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
3998   PlotSimple2D( canvas2DCorr, hspectraLGvsCellID, -10000, -10000, textSizeRel, Form("%s/LG.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
3999   PlotSimple2D( canvas2DCorr, hspectraHGCorrvsCellID, -10000, -10000, textSizeRel, Form("%s/HGCorr.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
4000   PlotSimple2D( canvas2DCorr, hspectraHGCorrvsCellID, 300, -10000, textSizeRel, Form("%s/HGCorr_zoomed.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
4001   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");
4002   PlotSimple2D( canvas2DCorr, hspectraLGCorrvsCellID, -10000, -10000, textSizeRel, Form("%s/LGCorr.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
4003   PlotSimple2D( canvas2DCorr, hspectraLGCorrvsCellID, 200, -10000, textSizeRel, Form("%s/LGCorr_zoomed.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
4004   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");
4005   PlotSimple2D( canvas2DCorr, hspectraEnergyvsCellID, -10000, -10000, textSizeRel, Form("%s/EnergyVsCellID.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
4006   PlotSimple2D( canvas2DCorr, hspectraEnergyTotvsNCells, -10000, -10000, textSizeRel, Form("%s/EnergyTotalVsNCells.%s", outputDirPlots.Data(), plotSuffix.Data()), it->second, 1, kFALSE, "colz", true);
4007   
4008   TCanvas* canvas1DSimple = new TCanvas("canvas1DSimple","",0,0,1450,1300);  // gives the page size
4009   DefaultCancasSettings( canvas1DSimple, 0.08, 0.03, 0.03, 0.07);
4010   hspectraEnergyTot->Scale(1./evts);
4011   hspectraEnergyTot->GetYaxis()->SetTitle("counts/event");
4012   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() ));
4013   hspectraNCells->Scale(1./evts);
4014   hspectraNCells->GetYaxis()->SetTitle("counts/event");
4015   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() ));
4016   
4017   if (ExtPlot > 0){
4018     //***********************************************************************************************************
4019     //********************************* 8 Panel overview plot  **************************************************
4020     //***********************************************************************************************************
4021     //*****************************************************************
4022       // Test beam geometry (beam coming from viewer)
4023       //===========================================================
4024       //||    8 (4)    ||    7 (5)   ||    6 (6)   ||    5 (7)   ||  row 0
4025       //===========================================================
4026       //||    1 (0)    ||    2 (1)   ||    3 (2)   ||    4 (3)   ||  row 1
4027       //===========================================================
4028       //    col 0     col 1       col 2     col  3
4029       // rebuild pad geom in similar way (numbering -1)
4030     //*****************************************************************
4031     TCanvas* canvas8Panel;
4032     TPad* pad8Panel[8];
4033     Double_t topRCornerX[8];
4034     Double_t topRCornerY[8];
4035     Int_t textSizePixel = 30;
4036     Double_t relSize8P[8];
4037     CreateCanvasAndPadsFor8PannelTBPlot(canvas8Panel, pad8Panel,  topRCornerX, topRCornerY, relSize8P, textSizePixel);
4038 
4039     TCanvas* canvas8PanelProf;
4040     TPad* pad8PanelProf[8];
4041     Double_t topRCornerXProf[8];
4042     Double_t topRCornerYProf[8];
4043     Double_t relSize8PProf[8];
4044     CreateCanvasAndPadsFor8PannelTBPlot(canvas8PanelProf, pad8PanelProf,  topRCornerXProf, topRCornerYProf, relSize8PProf, textSizePixel, 0.045, "Prof", false);
4045 
4046     calib.PrintGlobalInfo();  
4047     std::cout << "plotting single layers" << std::endl;
4048 
4049     for (Int_t l = 0; l < setup->GetNMaxLayer()+1; l++){
4050       for (Int_t m = 0; m < setup->GetNMaxModule()+1; m++){    
4051         if (l%10 == 0 && l > 0 && debug > 0)
4052           std::cout << "============================== layer " <<  l << " / " << setup->GetNMaxLayer() << " layers" << std::endl;
4053         if (ExtPlot > 0){
4054           PlotNoiseAdvWithFits8MLayer (canvas8Panel,pad8Panel, topRCornerX, topRCornerY, relSize8P, textSizePixel, 
4055                                       hSpectra, hSpectraNoise, true, -50, 100, 1.2, l, m,
4056                                       Form("%s/NoiseTrigg_HG_Mod%02d_Layer%02d.%s" ,outputDirPlots.Data(), m, l, plotSuffix.Data()), it->second);
4057           PlotNoiseAdvWithFits8MLayer (canvas8Panel,pad8Panel, topRCornerX, topRCornerY, relSize8P, textSizePixel, 
4058                                       hSpectra, hSpectraNoise, false, -50, 100, 1.2, l, m,
4059                                       Form("%s/NoiseTrigg_LG_Mod%02d_Layer%02d.%s" ,outputDirPlots.Data(), m, l, plotSuffix.Data()), it->second);
4060           PlotSpectra8MLayer (canvas8Panel,pad8Panel, topRCornerX, topRCornerY, relSize8P, textSizePixel, 
4061                                     hSpectra, 0, -100, 4000, 1.2, l, m,
4062                                     Form("%s/Spectra_HG_Mod%02d_Layer%02d.%s" ,outputDirPlots.Data(), m, l, plotSuffix.Data()), it->second);
4063           PlotSpectra8MLayer (canvas8Panel,pad8Panel, topRCornerX, topRCornerY, relSize8P, textSizePixel, 
4064                                     hSpectra, 2, -2, 100, 1.2, l, m,
4065                                     Form("%s/Spectra_Comb_Mod%02d_Layer%02d.%s" ,outputDirPlots.Data(), m, l, plotSuffix.Data()), it->second);
4066           PlotCorrWithFits8MLayer(canvas8PanelProf,pad8PanelProf, topRCornerXProf, topRCornerYProf, relSize8PProf, textSizePixel, 
4067                                     hSpectra, 0, -20, 800, 50., l, m,
4068                                     Form("%s/LGHG_Corr_Mod%02d_Layer%02d.%s" ,outputDirPlots.Data(), m, l, plotSuffix.Data()), it->second);
4069           PlotCorrWithFits8MLayer(canvas8PanelProf,pad8PanelProf, topRCornerXProf, topRCornerYProf, relSize8PProf, textSizePixel, 
4070                                     hSpectra, 2, -100, 340, 300., l, m,
4071                                     Form("%s/LGLGhgeq_Corr_Mod%02d_Layer%02d.%s" ,outputDirPlots.Data(), m, l, plotSuffix.Data()), it->second);
4072         }
4073       }
4074     }
4075     std::cout << "done plotting single layers" << std::endl;  
4076   }
4077   
4078   return true;
4079 }
4080 
4081 
4082 //***********************************************************************************************
4083 //*********************** Save Noise triggers only ***************************************************
4084 //***********************************************************************************************
4085 bool Analyses::SaveNoiseTriggersOnly(void){
4086   std::cout<<"Save noise triggers into separate file"<<std::endl;
4087   TcalibIn->GetEntry(0);
4088   // check whether calib should be overwritten based on external text file
4089   if (OverWriteCalib){
4090     calib.ReadCalibFromTextFile(ExternalCalibFile,debug);
4091   }
4092   
4093   int evts=TdataIn->GetEntries();
4094   for(int i=0; i<evts; i++){
4095     TdataIn->GetEntry(i);
4096     if (i%5000 == 0 && i > 0 && debug > 0) std::cout << "Reading " <<  i << " / " << evts << " events" << std::endl;
4097     for(int j=0; j<event.GetNTiles(); j++){
4098       Caen* aTile=(Caen*)event.GetTile(j);
4099       // testing for noise trigger
4100       if(aTile->GetLocalTriggerBit()!= (char)2){
4101         event.RemoveTile(aTile);
4102         j--;
4103       }
4104     }
4105     RootOutput->cd();
4106     TdataOut->Fill();
4107   }
4108   TdataOut->Write();
4109   TsetupIn->CloneTree()->Write();
4110   
4111   if (IsCalibSaveToFile()){
4112     TString fileCalibPrint = RootOutputName;
4113     fileCalibPrint         = fileCalibPrint.ReplaceAll(".root","_calib.txt");
4114     calib.PrintCalibToFile(fileCalibPrint);
4115   }
4116 
4117   TcalibOut->Fill();
4118   TcalibOut->Write();
4119   RootOutput->Close();
4120   RootInput->Close();      
4121   
4122   return true;
4123 }
4124 
4125 //***********************************************************************************************
4126 //*********************** Save Noise triggers only ***************************************************
4127 //***********************************************************************************************
4128 bool Analyses::SaveCalibToOutputOnly(void){
4129   std::cout<<"Save calib into separate file: "<< GetRootCalibOutputName().Data() <<std::endl;
4130   RootCalibOutput->cd();
4131   TcalibIn->GetEntry(0);  
4132   TsetupIn->CloneTree()->Write();
4133   
4134   // check whether calib should be overwritten based on external text file
4135   if (OverWriteCalib){
4136     calib.ReadCalibFromTextFile(ExternalCalibFile,debug);
4137   }
4138   
4139   if (IsCalibSaveToFile()){
4140     TString fileCalibPrint = GetRootCalibOutputName();
4141     fileCalibPrint         = fileCalibPrint.ReplaceAll(".root","_calib.txt");
4142     calib.PrintCalibToFile(fileCalibPrint);
4143   }
4144 
4145   TcalibOut->Fill();
4146   TcalibOut->Write();
4147   RootCalibOutput->Close();
4148   
4149   return true;
4150 }
4151 
4152 
4153 //***********************************************************************************************
4154 //*********************** Save local muon triggers only ***************************************************
4155 //***********************************************************************************************
4156 bool Analyses::SaveMuonTriggersOnly(void){
4157   std::cout<<"Save local muon triggers into separate file"<<std::endl;
4158   TcalibIn->GetEntry(0);
4159     // check whether calib should be overwritten based on external text file
4160   if (OverWriteCalib){
4161     calib.ReadCalibFromTextFile(ExternalCalibFile,debug);
4162   }
4163 
4164   int evts=TdataIn->GetEntries();
4165   for(int i=0; i<evts; i++){
4166     TdataIn->GetEntry(i);
4167     if (i%5000 == 0 && i > 0 && debug > 0) std::cout << "Reading " <<  i << " / " << evts << " events" << std::endl;
4168     for(int j=0; j<event.GetNTiles(); j++){
4169       Caen* aTile=(Caen*)event.GetTile(j);
4170       // testing for muon trigger
4171       if(aTile->GetLocalTriggerBit()!= (char)1){
4172         event.RemoveTile(aTile);
4173         j--;
4174       }
4175     }
4176     RootOutput->cd();
4177     TdataOut->Fill();
4178   }
4179   TdataOut->Write();
4180   TsetupIn->CloneTree()->Write();
4181   
4182   if (IsCalibSaveToFile()){
4183     TString fileCalibPrint = RootOutputName;
4184     fileCalibPrint         = fileCalibPrint.ReplaceAll(".root","_calib.txt");
4185     calib.PrintCalibToFile(fileCalibPrint);
4186   }
4187 
4188   TcalibOut->Fill();
4189   TcalibOut->Write();
4190   RootOutput->Close();
4191   RootInput->Close();      
4192   
4193   return true;
4194 }
4195 
4196 
4197 //***********************************************************************************************
4198 //*********************** Skim HGCROC data ******************************************************
4199 //***********************************************************************************************
4200 bool Analyses::SkimHGCROCData(void){
4201   std::cout<<"Skim HGCROC data from pure noise"<<std::endl;
4202   TcalibIn->GetEntry(0);
4203     // check whether calib should be overwritten based on external text file
4204   if (OverWriteCalib){
4205     calib.ReadCalibFromTextFile(ExternalCalibFile,debug);
4206   }
4207 
4208   int evts=TdataIn->GetEntries();
4209   if (maxEvents == -1){
4210     maxEvents = evts;
4211   } else {
4212     std::cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl;
4213     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;
4214     std::cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl;
4215   }
4216 
4217   long evtTrigg = 0;
4218   
4219   for(int i=0; i<evts && i < maxEvents; i++){
4220     TdataIn->GetEntry(i);
4221     if (i%5000 == 0 && i > 0 && debug > 0) std::cout << "Reading " <<  i << " / " << evts << " events" << std::endl;
4222     // std::cout << "Reading " <<  i << " / " << evts << " events" << std::endl;
4223     bool triggered = false;
4224     for(int j=0; j<event.GetNTiles(); j++){
4225       Hgcroc* aTile=(Hgcroc*)event.GetTile(j);
4226       // testing for any signal beyond noise
4227       // std::cout << "Cell ID:" << aTile->GetCellID() << "\t" << calib.GetPedestalMeanH(aTile->GetCellID()) << "\t" << calib.GetPedestalMeanL(aTile->GetCellID()) << "\n \tADC-wave " ;
4228       // for (int k = 0; k < (int)aTile->GetADCWaveform().size(); k++ ){
4229       //   std::cout << aTile->GetADCWaveform().at(k) << "\t" ;
4230       // }
4231       // // std::cout << "\n \tTOT-Wave ";
4232       // // for (int k = 0; k < (int)aTile->GetTOTWaveform().size(); k++ ){
4233       // //   std::cout << aTile->GetTOTWaveform().at(k) << "\t" ;
4234       // // }
4235       // std::cout << "\n \tTOA-Wave ";
4236       // for (int k = 0; k < (int)aTile->GetTOAWaveform().size(); k++ ){
4237       //   std::cout << aTile->GetTOAWaveform().at(k) << "\t" ;
4238       // }
4239       // std::cout <<"\n\t\t\t";
4240       // for (int k = 0; k < (int)aTile->GetTOAWaveform().size(); k++ )
4241       //   std::cout <<"\t";  
4242       // std::cout << " integ: "<< aTile->GetTOT() << "\t" << aTile->GetTOA() << std::endl;
4243       
4244       if (aTile->GetTOA() > 0) triggered= true;
4245       // if(aTile->GetLocalTriggerBit()!= (char)1){
4246         // event.RemoveTile(aTile);
4247         // j--;
4248       // }
4249     }
4250     
4251     if (!triggered && debug == 3){
4252       for(int j=0; j<event.GetNTiles(); j++){
4253         Hgcroc* aTile=(Hgcroc*)event.GetTile(j);
4254         // testing for any signal beyond noise
4255         std::cout << "Cell ID:" << aTile->GetCellID() << "\t" << calib.GetPedestalMeanH(aTile->GetCellID()) << "\t" << calib.GetPedestalMeanL(aTile->GetCellID()) << "\n \tADC-wave " ;
4256         for (int k = 0; k < (int)aTile->GetADCWaveform().size(); k++ ){
4257           std::cout << aTile->GetADCWaveform().at(k) << "\t" ;
4258         }
4259         std::cout << "\n \tTOT-Wave ";
4260         for (int k = 0; k < (int)aTile->GetTOTWaveform().size(); k++ ){
4261           std::cout << aTile->GetTOTWaveform().at(k) << "\t" ;
4262         }
4263         std::cout << "\n \tTOA-Wave ";
4264         for (int k = 0; k < (int)aTile->GetTOAWaveform().size(); k++ ){
4265           std::cout << aTile->GetTOAWaveform().at(k) << "\t" ;
4266         }
4267         std::cout <<"\n\t\t\t";
4268         for (int k = 0; k < (int)aTile->GetTOAWaveform().size(); k++ )
4269           std::cout <<"\t";  
4270         std::cout << " integ: "<< aTile->GetTOT() << "\t" << aTile->GetTOA() << std::endl;
4271       }
4272     }
4273     
4274     RootOutput->cd();
4275     if (triggered){
4276       evtTrigg++;
4277       TdataOut->Fill();
4278     }
4279   }
4280   TdataOut->Write();
4281   TsetupIn->CloneTree()->Write();
4282   
4283   std::cout << "Evts in: " << maxEvents << "\t skimmed: " << evtTrigg << std::endl;
4284    
4285   if (IsCalibSaveToFile()){
4286     TString fileCalibPrint = RootOutputName;
4287     fileCalibPrint         = fileCalibPrint.ReplaceAll(".root","_calib.txt");
4288     calib.PrintCalibToFile(fileCalibPrint);
4289   }
4290 
4291   TcalibOut->Fill();
4292   TcalibOut->Write();
4293   RootOutput->Close();
4294   RootInput->Close();      
4295   
4296   return true;
4297 }
4298 
4299 
4300 //***********************************************************************************************
4301 //*********************** Create output files ***************************************************
4302 //***********************************************************************************************
4303 bool Analyses::CreateOutputRootFile(void){
4304   if(Overwrite){
4305     std::cout << "overwriting exisiting output file" << std::endl;
4306     RootOutput=new TFile(RootOutputName.Data(),"RECREATE");
4307   } else{
4308     std::cout << "creating output file" << std::endl;
4309     RootOutput = new TFile(RootOutputName.Data(),"CREATE");
4310   }
4311   if(RootOutput->IsZombie()){
4312     std::cout<<"Error opening '"<<RootOutput<<"'no reachable path? Exist without force mode to overwrite?..."<<std::endl;
4313     return false;
4314   }
4315   return true;
4316 }
4317 
4318 //***********************************************************************************************
4319 //*********************** Read external bad channel map *****************************************
4320 //***********************************************************************************************
4321 std::map<int,short> Analyses::ReadExternalBadChannelMap(void){
4322   
4323   std::cout << "Reading in external mapping file" << std::endl;
4324   std::map<int,short> bcmap;
4325   
4326   std::ifstream bcmapFile;
4327   bcmapFile.open(ExternalBadChannelMap,std::ios_base::in);
4328   if (!bcmapFile) {
4329     std::cout << "ERROR: file " << ExternalBadChannelMap.Data() << " not found!" << std::endl;
4330     return bcmap;
4331   }
4332 
4333   for( TString tempLine; tempLine.ReadLine(bcmapFile, kTRUE); ) {
4334     // check if line should be considered
4335     if (tempLine.BeginsWith("%") || tempLine.BeginsWith("#")){
4336       continue;
4337     }
4338     if (debug > 1) std::cout << tempLine.Data() << std::endl;
4339 
4340     // Separate the string according to tabulators
4341     TObjArray *tempArr  = tempLine.Tokenize(" ");
4342     if(tempArr->GetEntries()<2){
4343       if (debug > 1) std::cout << "nothing to be done" << std::endl;
4344       delete tempArr;
4345       continue;
4346     } 
4347     
4348     int mod     = ((TString)((TObjString*)tempArr->At(0))->GetString()).Atoi();
4349     int layer   = ((TString)((TObjString*)tempArr->At(1))->GetString()).Atoi();
4350     int row     = ((TString)((TObjString*)tempArr->At(2))->GetString()).Atoi();
4351     int col     = ((TString)((TObjString*)tempArr->At(3))->GetString()).Atoi();
4352     short bc    = short(((TString)((TObjString*)tempArr->At(4))->GetString()).Atoi());
4353     
4354     int cellID  = setup->GetCellID( row, col, layer, mod);    
4355                 
4356     if (debug > 1) std::cout << "cellID " << cellID << "\t BC status: " << bc<< std::endl;
4357     bcmap[cellID]=bc;
4358   }
4359   std::cout << "registered " << bcmap.size() << " bad channels!" << std::endl;
4360   return bcmap;  
4361   
4362 }
[ 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!