epic-lfhcal-tbana/NewStructure/Analyses.cc

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