epic-lfhcal-tbana/NewStructure/Analyses.cc

0001 #include "Analyses.h"
0002 #include <vector>
0003 #include "TROOT.h"
0004 //#include <unistd.h> // Add for use on Mac OS -> Same goes for Analyse.cc
0005 #include "TF1.h"
0006 #include "TFitResult.h"
0007 #include "TFitResultPtr.h"
0008 #include "TH1D.h"
0009 #include "TH2D.h"
0010 #include "TProfile.h"
0011 #include "TChain.h"
0012 #include "CommonHelperFunctions.h"
0013 #include "TileSpectra.h"
0014 #include "PlottHelper.h"
0015 #include "TRandom3.h"
0016 #include "TMath.h"
0017 #include "Math/Minimizer.h"
0018 #include "Math/MinimizerOptions.h"
0019
0020 // ****************************************************************************
0021 // Checking and opening input and output files
0022 // ****************************************************************************
0023 bool Analyses::CheckAndOpenIO(void){
0024   int matchingbranch;
0025   if(!ASCIIinputName.IsNull()){
0026     ASCIIinput.open(ASCIIinputName.Data(),std::ios::in);
0027     if(!ASCIIinput.is_open()){
0028       std::cout<<"Could not open input file: "<<optarg<<std::endl;
0029       return false;
0030     }
0031   }
0032
0033   //Need to check first input to get the setup...I do not think it is necessary
0034     std::cout<<"Input name set to: '"<<RootInputName.Data() <<std::endl;
0035
0036   if(!RootInputName.IsNull()){
0037     //File exist?
0038     RootInput=new TFile(RootInputName.Data(),"READ");
0039     if(RootInput->IsZombie()){
0040       std::cout<<"Error opening '"<<RootInputName<<"', does the file exist?"<<std::endl;
0041       return false;
0042     }
0043
0044     //Retrieve info, start with setup
0045     TsetupIn = (TTree*) RootInput->Get("Setup");
0046     if(!TsetupIn){
0047       std::cout<<"Could not retrieve the Setup tree, leaving"<<std::endl;
0048       return false;
0049     }
0050     setup=Setup::GetInstance();
0051     std::cout<<"Setup add "<<setup<<std::endl;
0052     //matchingbranch=TsetupIn->SetBranchAddress("setup",&setup);
0053     matchingbranch=TsetupIn->SetBranchAddress("setup",&rswptr);
0054     if(matchingbranch<0){
0055       std::cout<<"Error retrieving Setup info from the tree"<<std::endl;
0056       return false;
0057     }
0058     std::cout<<"Entries "<<TsetupIn->GetEntries()<<std::endl;
0059     TsetupIn->GetEntry(0);
0060     setup->Initialize(*rswptr);
0061     std::cout<<"Reading "<<RootInput->GetName()<<std::endl;
0062     std::cout<<"Setup Info "<<setup->IsInit()<<"  and  "<<setup->GetCellID(0,0)<<std::endl;
0063     //std::cout<<"Setup add now "<<setup<<std::endl;
0064
0065     //Retrieve info, existing data?
0066     TdataIn = (TTree*) RootInput->Get("Data");
0067     if(!TdataIn){
0068       std::cout<<"Could not retrieve the Data tree, leaving"<<std::endl;
0069       return false;
0070     }
0071     matchingbranch=TdataIn->SetBranchAddress("event",&eventptr);
0072     if(matchingbranch<0){
0073       std::cout<<"Error retrieving Event info from the tree"<<std::endl;
0074       return false;
0075     }
0076
0077     //Do I really want this?
0078     TcalibIn = (TTree*) RootInput->Get("Calib");
0079     if(!TcalibIn){
0080       std::cout<<"Could not retrieve Calib tree, leaving"<<std::endl;
0081       //return false;
0082     }
0083     else {
0084       matchingbranch=TcalibIn->SetBranchAddress("calib",&calibptr);
0085       if(matchingbranch<0){
0086         std::cout<<"Error retrieving calibration info from the tree"<<std::endl;
0087         TcalibIn=nullptr;
0088       }
0089     }
0090     //End of do I really want this?
0091
0092     //We want to retrieve also calibration if do not specify ApplyPedestalCorrection from external file
0093     //In other words, the pedestal was potentially already done and we have an existing calib object
0094     if((!ApplyPedestalCorrection && ExtractScaling) || (!ApplyPedestalCorrection && ExtractScalingImproved) || (!ApplyPedestalCorrection && ReextractNoise)){
0095       TcalibIn = (TTree*) RootInput->Get("Calib");
0096       if(!TcalibIn){
0097         std::cout<<"Could not retrieve Calib tree, leaving"<<std::endl;
0098         return false;
0099       }
0100       //calib=Calib::GetInstance();
0101       matchingbranch=TcalibIn->SetBranchAddress("calib",&calibptr);
0102       if(matchingbranch<0){
0103         std::cout<<"Error retrieving calibration info from the tree"<<std::endl;
0104         return false;
0105       }
0106     }
0107     //All good
0108   }
0109   else if(!Convert){
0110     std::cout<<"Explicit Input option mandatory except for convertion ASCII -> ROOT"<<std::endl;
0111     return false;
0112   }
0113
0114   //Setup Output files
0115   if(!RootOutputName.IsNull()){
0116     if (!Convert){
0117       TString temp = RootOutputName;
0118       temp         = temp.ReplaceAll(".root","_Hists.root");
0119       SetRootOutputHists(temp);
0120       std::cout << "creating additional histo file: " << RootOutputNameHist.Data() << " tree in : "<< RootOutputName.Data() << std::endl;
0121     }
0122
0123     bool sCOF = CreateOutputRootFile();
0124     if (!sCOF) return false;
0125
0126     TsetupOut = new TTree("Setup","Setup");
0127     setup=Setup::GetInstance();
0128     //TsetupOut->Branch("setup",&setup);
0129     TsetupOut->Branch("setup",&rsw);
0130     TdataOut = new TTree("Data","Data");
0131     TdataOut->Branch("event",&event);
0132     //if(!calib) calib=new Calib();
0133     //if(!calib)
0134     //calib=Calib::GetInstance();
0135     //Calib* calib=Calib::GetInstance();
0136     //std::cout<<"Calib pointer is "<<calibptr<<std::endl;
0137     TcalibOut = new TTree("Calib","Calib");
0138     TcalibOut->Branch("calib",&calib);
0139   }
0140   else {
0141     std::cout<<"Output option mandatory except when converting"<<std::endl;
0142     return false;
0143   }
0144
0145   if(!RootCalibInputName.IsNull()){
0146     RootCalibInput=new TFile(RootCalibInputName.Data(),"READ");
0147     if(RootCalibInput->IsZombie()){
0148       std::cout<<"Error opening '"<<RootCalibInputName<<"', does the file exist?"<<std::endl;
0149       return false;
0150     }
0151     TcalibIn = nullptr;
0152     TcalibIn = (TTree*) RootCalibInput->Get("Calib");
0153     if(!TcalibIn){
0154       std::cout<<"Could not retrieve Calib tree, leaving"<<std::endl;
0155       return false;
0156     } else {
0157       std::cout<<"Retrieved calib object from external input!"<<std::endl;
0158     }
0159     matchingbranch=TcalibIn->SetBranchAddress("calib",&calibptr);
0160     if(matchingbranch<0){
0161       std::cout<<"Error retrieving calibration info from the tree"<<std::endl;
0162       return false;
0163     }else {
0164       std::cout<<"Correctly set branch for external calib input."<<std::endl;
0165     }
0166
0167   }
0168
0169   if(!RootPedestalInputName.IsNull()){
0170     RootPedestalInput = new TFile(RootPedestalInputName.Data(),"READ");
0171     if(RootPedestalInput->IsZombie()){
0172       std::cout<<"Error opening '"<<RootPedestalInputName<<"', does the file exist?"<<std::endl;
0173       return false;
0174     }
0175     TcalibIn = (TTree*) RootPedestalInput->Get("Calib");
0176     if(!TcalibIn){
0177       std::cout<<"Could not retrieve Calib tree, leaving"<<std::endl;
0178       return false;
0179     } else {
0180       std::cout<<"Retrieved calib object from external input for pedestals!"<<std::endl;
0181     }
0182     matchingbranch=TcalibIn->SetBranchAddress("calib",&calibptr);
0183     if(matchingbranch<0){
0184       std::cout<<"Error retrieving calibration info from the tree"<<std::endl;
0185       return false;
0186     }else {
0187       std::cout<<"Correctly set branch for external calib for pedestal input."<<std::endl;
0188     }
0189     //std::cout<<"Did the address changed? "<<&calib<<std::endl;
0190   }
0191   if(!MapInputName.IsNull()){
0192     MapInput.open(MapInputName.Data(),std::ios::in);
0193     if(!MapInput.is_open()){
0194       std::cout<<"Could not open mapping file: "<<MapInputName<<std::endl;
0195       return false;
0196     }
0197   }
0198   return true;
0199 }
0200
0201 // ****************************************************************************
0202 // Primary process function to start all calibrations
0203 // ****************************************************************************
0204 bool Analyses::Process(void){
0205   bool status;
0206   ROOT::EnableImplicitMT();
0207
0208   if(Convert){
0209     if (!(GetASCIIinputName().EndsWith(".root"))){
0210       status=ConvertASCII2Root();
0211     } else {
0212       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;
0213       status=ConvertOldRootFile2Root();
0214     }
0215   }
0216
0217   // extract pedestal from pure pedestal runs (internal triggers)
0218   if(ExtractPedestal){
0219     status=GetPedestal();
0220   }
0221
0222   // copy existing pedestal to other file
0223   if(ApplyPedestalCorrection){
0224     status=CorrectPedestal();
0225   }
0226
0227   // extract mip calibration
0228   if(ExtractScaling){
0229     status=GetScaling();
0230   }
0231
0232   // extract noise sample from trigger flagged files
0233   if(ReextractNoise){
0234     status=GetNoiseSampleAndRefitPedestal();
0235   }
0236
0237   // rerun mip calibration based on triggers
0238   if (ExtractScalingImproved){
0239     status=GetImprovedScaling();
0240   }
0241
0242   // copy full calibration to different file and calculate energy
0243   if(ApplyCalibration){
0244     status=Calibrate();
0245   }
0246
0247   // reduce file to only noise triggers
0248   if(SaveNoiseOnly){
0249     status=SaveNoiseTriggersOnly();
0250   }
0251
0252   // reduce file to only mip triggers
0253   if(SaveMipsOnly){
0254     status=SaveMuonTriggersOnly();
0255   }
0256
0257   return status;
0258 }
0259
0260
0261 // ****************************************************************************
0262 // convert original ASCII file from CAEN output to root format
0263 // ****************************************************************************
0264 bool Analyses::ConvertASCII2Root(void){
0265   //============================================
0266   //Init first
0267   //============================================
0268   // initialize setup
0269   if (MapInputName.CompareTo("")== 0) {
0270       std::cerr << "ERROR: No mapping file has been provided, please make sure you do so! Aborting!" << std::endl;
0271       return false;
0272   }
0273   setup->Initialize(MapInputName.Data(),debug);
0274   // initialize run number file
0275   if (RunListInputName.CompareTo("")== 0) {
0276       std::cerr << "ERROR: No run list file has been provided, please make sure you do so! Aborting!" << std::endl;
0277       return false;
0278   }
0279   std::map<int,RunInfo> ri=readRunInfosFromFile(RunListInputName.Data(),debug,0);
0280
0281   // Clean up file headers
0282   TObjArray* tok=ASCIIinputName.Tokenize("/");
0283   // get file name
0284   TString file=((TObjString*)tok->At(tok->GetEntries()-1))->String();
0285   delete tok;
0286   tok=file.Tokenize("_");
0287   TString header=((TObjString*)tok->At(0))->String();
0288   delete tok;
0289
0290   // Get run number from file & find necessary run infos
0291   TString RunString=header(3,header.Sizeof());
0292   std::map<int,RunInfo>::iterator it=ri.find(RunString.Atoi());
0293   //std::cout<<RunString.Atoi()<<"\t"<<it->second.species<<std::endl;
0294   event.SetRunNumber(RunString.Atoi());
0295   event.SetROtype(2);
0296   event.SetBeamName(it->second.species);
0297   event.SetBeamID(it->second.pdg);
0298   event.SetBeamEnergy(it->second.energy);
0299   event.SetVop(it->second.vop);
0300   event.SetVov(it->second.vop-it->second.vbr);
0301   event.SetBeamPosX(it->second.posX);
0302   event.SetBeamPosY(it->second.posY);
0303   calib.SetRunNumber(RunString.Atoi());
0304   calib.SetVop(it->second.vop);
0305   calib.SetVov(it->second.vop-it->second.vbr);
0306
0307   //============================================
0308   // Start decoding file
0309   //============================================
0310   TString aline         = "";
0311   TString versionCAEN   = "";
0312   TObjArray* tokens;
0313   std::map<int,std::vector<Caen> > tmpEvent;
0314   std::map<int,double> tmpTime;
0315   std::map<int,std::vector<Caen> >::iterator itevent;
0316   long tempEvtCounter = 0;
0317   while(!ASCIIinput.eof()){                                                     // run till end of file is reached and read line by line
0318     aline.ReadLine(ASCIIinput);
0319     if(!ASCIIinput.good()) break;
0320     if(aline[0]=='/'){
0321       if (aline.Contains("File Format Version")){
0322         tokens      = aline.Tokenize(" ");
0323         versionCAEN = ((TObjString*)tokens->At(4))->String();
0324         std::cout << "File version: " << ((TObjString*)tokens->At(4))->String() << std::endl;
0325
0326         if (!(versionCAEN.CompareTo("3.3") == 0 || versionCAEN.CompareTo("3.1") == 0)){
0327           std::cerr << "This version cannot be converted with the current code, please add the corresponding file format to the converter" << std::endl;
0328           tokens->Clear();
0329           delete tokens;
0330           return false;
0331         }
0332
0333         tokens->Clear();
0334         delete tokens;
0335       }
0336       else if(aline.Contains("Run start time")){
0337         tokens    = aline.Tokenize(" ");
0338         int year=((TObjString*)tokens->At(8))->String().Atoi();
0339         int month;
0340         TString Stringmonth=((TObjString*)tokens->At(5))->String();
0341         if(Stringmonth=="Jan") month=1;
0342         else if(Stringmonth=="Feb") month=2;
0343         else if(Stringmonth=="Mar") month=3;
0344         else if(Stringmonth=="Apr") month=4;
0345         else if(Stringmonth=="May") month=5;
0346         else if(Stringmonth=="Jun") month=6;
0347         else if(Stringmonth=="Jul") month=7;
0348         else if(Stringmonth=="Aug") month=8;
0349         else if(Stringmonth=="Sep") month=9;
0350         else if(Stringmonth=="Oct") month=10;
0351         else if(Stringmonth=="Nov") month=11;
0352         else if(Stringmonth=="Dec") month=12;
0353         int day=((TObjString*)tokens->At(6))->String().Atoi();
0354         int hour=((TString)((TObjString*)tokens->At(7))->String()(0,2)).Atoi();
0355         int min=((TString)((TObjString*)tokens->At(7))->String()(3,5)).Atoi();
0356         int sec=((TString)((TObjString*)tokens->At(7))->String()(6,8)).Atoi();
0357         TTimeStamp t=TTimeStamp(year,month,day,hour,min,sec);
0358         event.SetBeginRunTime(t);
0359         calib.SetBeginRunTime(t);
0360         tokens->Clear();
0361         delete tokens;
0362       }
0363       continue;
0364     }
0365     tokens=aline.Tokenize(" \t");
0366     tokens->SetOwner(true);
0367
0368     if (versionCAEN.CompareTo("3.3") == 0){
0369       int Nfields=tokens->GetEntries();
0370       if(((TObjString*)tokens->At(0))->String()=="Brd") {
0371         tokens->Clear();
0372         delete tokens;
0373         continue;
0374       }
0375       //================================================================================
0376       // data format example
0377       // Brd  Ch       LG       HG        Tstamp_us        TrgID        NHits
0378       // 7  00       51       68     10203578.136            0      64
0379       // 7  01       55       75
0380       //================================================================================
0381       if(Nfields!=7){
0382         std::cout<<"Unexpected number of fields"<<std::endl;
0383         std::cout<<"Expected 7, got: "<<Nfields<<", exit"<<std::endl;
0384         for(int j=0; j<Nfields;j++) std::cout<<j<<"  "<<((TObjString*)tokens->At(j))->String()<<std::endl;
0385         tokens->Clear();
0386         delete tokens;
0387         return -1;
0388       }
0389       int TriggerID=((TObjString*)tokens->At(5))->String().Atoi();
0390       int NHits=((TObjString*)tokens->At(6))->String().Atoi();
0391       double Time = ((TObjString*)tokens->At(4))->String().Atof();
0392       Caen aTile;
0393       int aBoard=((TObjString*)tokens->At(0))->String().Atoi();
0394       int achannel=((TObjString*)tokens->At(1))->String().Atoi();
0395       aTile.SetE(((TObjString*)tokens->At(3))->String().Atoi());//To Test
0396       aTile.SetADCHigh(((TObjString*)tokens->At(3))->String().Atoi());
0397       aTile.SetADCLow (((TObjString*)tokens->At(2))->String().Atoi());
0398       tokens->Clear();
0399       delete tokens;
0400       aTile.SetCellID(setup->GetCellID(aBoard,achannel));
0401       itevent=tmpEvent.find(TriggerID);
0402
0403       if(itevent!=tmpEvent.end()){
0404         tmpTime[TriggerID]+=Time;
0405         if (aTile.GetCellID() != -1){
0406           itevent->second.push_back(aTile);
0407         } else {
0408           if(debug ==10) std::cout << "cell " << aBoard << "\t" << achannel << " wasn't active according to mapping file!" << std::endl;
0409         }
0410         for(int ich=1; ich<NHits; ich++){
0411             aline.ReadLine(ASCIIinput);
0412             tokens=aline.Tokenize(" ");
0413             tokens->SetOwner(true);
0414             Nfields=tokens->GetEntries();
0415             if(Nfields!=4){
0416               std::cout<<"Expecting 4 fields but read "<<Nfields<<std::endl;
0417               return -1;
0418             }
0419             achannel=((TObjString*)tokens->At(1))->String().Atoi();
0420             aTile.SetE(((TObjString*)tokens->At(3))->String().Atoi());//To Test
0421             aTile.SetADCHigh(((TObjString*)tokens->At(3))->String().Atoi());
0422             aTile.SetADCLow (((TObjString*)tokens->At(2))->String().Atoi());
0423             aTile.SetCellID(setup->GetCellID(aBoard,achannel));
0424             if (aTile.GetCellID() != -1){
0425               itevent->second.push_back(aTile);
0426             } else {
0427               if(debug ==10) std::cout << "cell " << aBoard << "\t" << achannel << " wasn't active according to mapping file!" << std::endl;
0428             }
0429             tokens->Clear();
0430             delete tokens;
0431         }
0432         if((int)itevent->second.size()==setup->GetTotalNbChannels()/*8*64*/){
0433           //Fill the tree the event is complete and erase from the map
0434           event.SetTimeStamp(tmpTime[TriggerID]/setup->GetNMaxROUnit());
0435           event.SetEventID(itevent->first);
0436           for(std::vector<Caen>::iterator itv=itevent->second.begin(); itv!=itevent->second.end(); ++itv){
0437             event.AddTile(new Caen(*itv));
0438           }
0439           TdataOut->Fill();
0440           tmpEvent.erase(itevent);
0441           tmpTime.erase(TriggerID);
0442         }
0443       } else {
0444         //This is a new event;
0445         tempEvtCounter++;                                                                   // in crease event counts for monitoring of progress
0446         if (tempEvtCounter%5000 == 0 && debug > 0) std::cout << "Converted " <<  tempEvtCounter << " events" << std::endl;
0447         std::vector<Caen> vCaen;
0448         if (aTile.GetCellID() != -1){
0449           vCaen.push_back(aTile);
0450         } else {
0451           if(debug ==10) std::cout << "cell " << aBoard << "\t" << achannel << " wasn't active according to mapping file!" << std::endl;
0452         }
0453         for(int ich=1; ich<NHits; ich++){
0454           aline.ReadLine(ASCIIinput);
0455           tokens=aline.Tokenize(" ");
0456           tokens->SetOwner(true);
0457           Nfields=tokens->GetEntries();
0458           if(Nfields!=4){
0459             std::cout<<"Expecting 4 fields but read "<<Nfields<<std::endl;
0460             return -1;
0461           }
0462           achannel=((TObjString*)tokens->At(1))->String().Atoi();
0463           aTile.SetE(((TObjString*)tokens->At(3))->String().Atoi());//To Test
0464           aTile.SetADCHigh(((TObjString*)tokens->At(3))->String().Atoi());
0465           aTile.SetADCLow (((TObjString*)tokens->At(2))->String().Atoi());
0466           aTile.SetCellID(setup->GetCellID(aBoard,achannel));
0467           if (aTile.GetCellID() != -1){
0468             vCaen.push_back(aTile);
0469           } else {
0470             if(debug ==10) std::cout << "cell " << aBoard << "\t" << achannel << " wasn't active according to mapping file!" << std::endl;
0471           }
0472           tokens->Clear();
0473           delete tokens;
0474         }
0475         tmpTime[TriggerID]=Time;
0476         tmpEvent[TriggerID]=vCaen;
0477       }
0478     } else if (versionCAEN.CompareTo("3.1") == 0){
0479       int Nfields=tokens->GetEntries();
0480       if(((TObjString*)tokens->At(0))->String()=="Tstamp_us") {
0481         tokens->Clear();
0482         delete tokens;
0483         continue;
0484       }
0485
0486       //================================================================================
0487       // data format example
0488       //   Tstamp_us        TrgID   Brd  Ch       LG       HG
0489       // 4970514.360            0    01  00       49       46
0490       //                             01  01       49       35
0491       //================================================================================
0492
0493       if(Nfields!=6){
0494         std::cout<<"Unexpected number of fields"<<std::endl;
0495         std::cout<<"Expected 6, got: "<<Nfields<<", exit"<<std::endl;
0496         for(int j=0; j<Nfields;j++) std::cout<<j<<"  "<<((TObjString*)tokens->At(j))->String()<<std::endl;
0497         tokens->Clear();
0498         delete tokens;
0499         return -1;
0500       }
0501
0502       // std::cout << aline.Data() << std::endl;
0503       int TriggerID = ((TObjString*)tokens->At(1))->String().Atoi();
0504       int NHits     = 64;                                           // temporary fix
0505       double Time   = ((TObjString*)tokens->At(0))->String().Atof();
0506       Caen aTile;
0507       int aBoard    =((TObjString*)tokens->At(2))->String().Atoi();
0508       int achannel  =((TObjString*)tokens->At(3))->String().Atoi();
0509       aTile.SetE(((TObjString*)tokens->At(5))->String().Atoi());//To Test
0510       aTile.SetADCHigh(((TObjString*)tokens->At(5))->String().Atoi());
0511       aTile.SetADCLow (((TObjString*)tokens->At(4))->String().Atoi());
0512       tokens->Clear();
0513       delete tokens;
0514       aTile.SetCellID(setup->GetCellID(aBoard,achannel));
0515       itevent=tmpEvent.find(TriggerID);
0516
0517       if(itevent!=tmpEvent.end()){
0518         tmpTime[TriggerID]+=Time;
0519         if (aTile.GetCellID() != -1){
0520           itevent->second.push_back(aTile);
0521         } else {
0522           if(debug ==10) std::cout << "cell " << aBoard << "\t" << achannel << " wasn't active according to mapping file!" << std::endl;
0523         }
0524         for(int ich=1; ich<NHits; ich++){
0525             aline.ReadLine(ASCIIinput);
0526             // std::cout << aline.Data() << std::endl;
0527             tokens=aline.Tokenize(" ");
0528             tokens->SetOwner(true);
0529             Nfields=tokens->GetEntries();
0530
0531             if(Nfields!=4){
0532               std::cout<< "Current line :" << aline.Data() << std::endl;
0533               std::cout<<"Expecting 4 fields but read "<<Nfields<<std::endl;
0534               return -1;
0535             }
0536             achannel=((TObjString*)tokens->At(1))->String().Atoi();
0537             aTile.SetE(((TObjString*)tokens->At(3))->String().Atoi());//To Test
0538             aTile.SetADCHigh(((TObjString*)tokens->At(3))->String().Atoi());
0539             aTile.SetADCLow (((TObjString*)tokens->At(2))->String().Atoi());
0540             aTile.SetCellID(setup->GetCellID(aBoard,achannel));
0541
0542             if (aTile.GetCellID() != -1){
0543               itevent->second.push_back(aTile);
0544             } else {
0545               if(debug ==10) std::cout << "cell " << aBoard << "\t" << achannel << " wasn't active according to mapping file!" << std::endl;
0546             }
0547             tokens->Clear();
0548             delete tokens;
0549         }
0550         if((int)itevent->second.size()==setup->GetTotalNbChannels()){
0551           //Fill the tree the event is complete and erase from the map
0552           event.SetTimeStamp(tmpTime[TriggerID]/setup->GetNMaxROUnit());
0553           event.SetEventID(itevent->first);
0554           for(std::vector<Caen>::iterator itv=itevent->second.begin(); itv!=itevent->second.end(); ++itv){
0555             event.AddTile(new Caen(*itv));
0556           }
0557           TdataOut->Fill();
0558           tmpEvent.erase(itevent);
0559           tmpTime.erase(TriggerID);
0560         }
0561       } else {
0562         //This is a new event;
0563         tempEvtCounter++;                                                                   // in crease event counts for monitoring of progress
0564         if (tempEvtCounter%5000 == 0 && debug > 0) std::cout << "Converted " <<  tempEvtCounter << " events" << std::endl;
0565         std::vector<Caen> vCaen;
0566
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         for(int ich=1; ich<NHits; ich++){
0573           aline.ReadLine(ASCIIinput);
0574           // std::cout << aline.Data() << std::endl;
0575           tokens=aline.Tokenize(" ");
0576           tokens->SetOwner(true);
0577           Nfields=tokens->GetEntries();
0578           if(Nfields!=4){
0579             std::cout<< "Current line :" << aline.Data() << std::endl;
0580             std::cout<<"Expecting 4 fields but read "<<Nfields<<std::endl;
0581             return -1;
0582           }
0583           achannel=((TObjString*)tokens->At(1))->String().Atoi();
0584           aTile.SetE(((TObjString*)tokens->At(3))->String().Atoi());//To Test
0585           aTile.SetADCHigh(((TObjString*)tokens->At(3))->String().Atoi());
0586           aTile.SetADCLow (((TObjString*)tokens->At(2))->String().Atoi());
0587           aTile.SetCellID(setup->GetCellID(aBoard,achannel));
0588           if (aTile.GetCellID() != -1){
0589             vCaen.push_back(aTile);
0590           } else {
0591             if(debug ==10) std::cout << "cell " << aBoard << "\t" << achannel << " wasn't active according to mapping file!" << std::endl;
0592           }
0593           tokens->Clear();
0594           delete tokens;
0595         }
0596         tmpTime[TriggerID]=Time;
0597         tmpEvent[TriggerID]=vCaen;
0598       }
0599     }
0600   } // finished reading in file
0601   //
0602   if (debug > 0) std::cout << "Converted a total of " <<  tempEvtCounter << " events" << std::endl;
0603
0604   //============================================
0605   // Fill & write all trees to output file
0606   //============================================
0607   RootOutput->cd();
0608   // setup
0609   RootSetupWrapper rswtmp=RootSetupWrapper(setup);
0610   rsw=rswtmp;
0611   TsetupOut->Fill();
0612   // calib
0613   TcalibOut->Fill();
0614   TcalibOut->Write();
0615   // event data
0616   TdataOut->Fill();
0617   TsetupOut->Write();
0618   TdataOut->Write();
0619
0620   RootOutput->Close();
0621   return true;
0622 }
0623
0624
0625
0626 // ****************************************************************************
0627 // convert already processed root file from CAEN output to new root format
0628 // ****************************************************************************
0629 bool Analyses::ConvertOldRootFile2Root(void){
0630   //============================================
0631   //Init first
0632   //============================================
0633   // initialize setup
0634   if (MapInputName.CompareTo("")== 0) {
0635       std::cerr << "ERROR: No mapping file has been provided, please make sure you do so! Aborting!" << std::endl;
0636       return false;
0637   }
0638   setup->Initialize(MapInputName.Data(),debug);
0639   // initialize run number file
0640   if (RunListInputName.CompareTo("")== 0) {
0641       std::cerr << "ERROR: No run list file has been provided, please make sure you do so! Aborting!" << std::endl;
0642       return false;
0643   }
0644   std::map<int,RunInfo> ri=readRunInfosFromFile(RunListInputName.Data(),debug,0);
0645
0646   // Clean up file headers
0647   TObjArray* tok=ASCIIinputName.Tokenize("/");
0648   // get file name
0649   TString file=((TObjString*)tok->At(tok->GetEntries()-1))->String();
0650   delete tok;
0651   tok=file.Tokenize("_");
0652   TString header=((TObjString*)tok->At(0))->String();
0653   delete tok;
0654
0655   // Get run number from file & find necessary run infos
0656   TString RunString=header(3,header.Sizeof());
0657   std::map<int,RunInfo>::iterator it=ri.find(RunString.Atoi());
0658   //std::cout<<RunString.Atoi()<<"\t"<<it->second.species<<std::endl;
0659   event.SetRunNumber(RunString.Atoi());
0660   event.SetROtype(2);
0661   event.SetBeamName(it->second.species);
0662   event.SetBeamID(it->second.pdg);
0663   event.SetBeamEnergy(it->second.energy);
0664   event.SetVop(it->second.vop);
0665   event.SetVov(it->second.vop-it->second.vbr);
0666   event.SetBeamPosX(it->second.posX);
0667   event.SetBeamPosY(it->second.posY);
0668   calib.SetRunNumber(RunString.Atoi());
0669   calib.SetVop(it->second.vop);
0670   calib.SetVov(it->second.vop-it->second.vbr);
0671
0672     // load tree
0673   TChain *const tt_event = new TChain("tree");
0674   if (ASCIIinputName.EndsWith(".root")) {                     // are we loading a single root tree?
0675       std::cout << "loading a single root file" << std::endl;
0676       tt_event->AddFile(ASCIIinputName);
0677       TFile testFile(ASCIIinputName.Data());
0678       if (testFile.IsZombie()){
0679         std::cout << Form("The file %s is not a root file or doesn't exit, please fix the file path", ASCIIinputName.Data()) << std::endl;
0680         return false;
0681       }
0682
0683   } else {
0684       std::cout << "please try again this isn't a root file" << std::endl;
0685   }
0686   if(!tt_event){ std::cout << "tree not found... returning!"<< std::endl; return false;}
0687
0688   // Define tree variables
0689   Long64_t gTrID;
0690   Double_t gTRtimeStamp;
0691   const int gMaxChannels = 64;
0692   Long64_t* gBoard          = new Long64_t[gMaxChannels];
0693   Long64_t* gChannel        = new Long64_t[gMaxChannels];
0694   Long64_t* gLG             = new Long64_t[gMaxChannels];
0695   Long64_t* gHG             = new Long64_t[gMaxChannels];
0696
0697   if (tt_event->GetBranchStatus("t_stamp") ){
0698     tt_event->SetBranchAddress("trgid",            &gTrID);
0699     tt_event->SetBranchAddress("t_stamp",          &gTRtimeStamp);
0700     tt_event->SetBranchAddress("board",            gBoard);
0701     tt_event->SetBranchAddress("channel",          gChannel);
0702     tt_event->SetBranchAddress("LG",               gLG);
0703     tt_event->SetBranchAddress("HG",               gHG);
0704   }
0705
0706   Long64_t nEntriesTree                 = tt_event->GetEntries();
0707   std::cout << "Number of events in tree: " << nEntriesTree << std::endl;
0708
0709   std::map<int,std::vector<Caen> > tmpEvent;
0710   std::map<int,double> tmpTime;
0711   for (Long64_t i=0; i<nEntriesTree;i++) {
0712     // load current event
0713     tt_event->GetEntry(i);
0714     if (i%5000 == 0 && debug > 0) std::cout << "Converted " <<  i << " events" << std::endl;
0715     int TriggerID = gTrID;
0716     double Time   = gTRtimeStamp;
0717     std::vector<Caen> vCaen;
0718
0719     for(Int_t ch=0; ch<gMaxChannels; ch++){
0720       Caen aTile;
0721       int aBoard      = gBoard[ch];
0722       int achannel    = gChannel[ch];
0723       aTile.SetE(gHG[ch]);//To Test
0724       aTile.SetADCHigh(gHG[ch]);
0725       aTile.SetADCLow (gLG[ch]);
0726       aTile.SetCellID(setup->GetCellID(aBoard,achannel));
0727       if (aTile.GetCellID() != -1){
0728         vCaen.push_back(aTile);
0729       } else {
0730         if(debug ==10) std::cout << "cell " << aBoard << "\t" << achannel << " wasn't active according to mapping file!" << std::endl;
0731       }
0732     }
0733
0734      if((int)vCaen.size()==setup->GetTotalNbChannels()){
0735       //Fill the tree the event is complete and erase from the map
0736       event.SetTimeStamp(Time);
0737       event.SetEventID(TriggerID);
0738       for(std::vector<Caen>::iterator itv=vCaen.begin(); itv!=vCaen.end(); ++itv){
0739         event.AddTile(new Caen(*itv));
0740       }
0741       TdataOut->Fill();
0742       vCaen.clear();
0743     }
0744   } // finished reading the file
0745
0746   //
0747   if (debug > 0) std::cout << "Converted a total of " <<  nEntriesTree << " events" << std::endl;
0748
0749   //============================================
0750   // Fill & write all trees to output file
0751   //============================================
0752   RootOutput->cd();
0753   // setup
0754   RootSetupWrapper rswtmp=RootSetupWrapper(setup);
0755   rsw=rswtmp;
0756   TsetupOut->Fill();
0757   // calib
0758   TcalibOut->Fill();
0759   TcalibOut->Write();
0760   // event data
0761   TdataOut->Fill();
0762   TsetupOut->Write();
0763   TdataOut->Write();
0764
0765   RootOutput->Close();
0766   return true;
0767 }
0768
0769
0770
0771 // ****************************************************************************
0772 // extract pedestral from dedicated data run, no other data in run
0773 // ****************************************************************************
0774 bool Analyses::GetPedestal(void){
0775   std::cout<<"GetPedestal for maximium "<< setup->GetMaxCellID() << " cells" <<std::endl;
0776
0777   std::map<int,RunInfo> ri=readRunInfosFromFile(RunListInputName.Data(),debug,0);
0778
0779   int maxChannelPerLayer = (setup->GetNMaxColumn()+1)*(setup->GetNMaxRow()+1);
0780   ROOT::Math::MinimizerOptions::SetDefaultMinimizer("Minuit2", "Migrad");
0781
0782   // create HG and LG histo's per channel
0783   TH2D* hspectraHGvsCellID      = new TH2D( "hspectraHG_vsCellID","ADC spectrum High Gain vs CellID; cell ID; ADC_{HG} (arb. units) ",
0784                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5, 4000,0,4000);
0785   hspectraHGvsCellID->SetDirectory(0);
0786   TH2D* hspectraLGvsCellID      = new TH2D( "hspectraLG_vsCellID","ADC spectrum Low Gain vs CellID; cell ID; ADC_{LG} (arb. units)  ",
0787                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5, 4000,0,4000);
0788   hspectraLGvsCellID->SetDirectory(0);
0789   TH1D* hMeanPedHGvsCellID      = new TH1D( "hMeanPedHG_vsCellID","mean Ped High Gain vs CellID ; cell ID; #mu_{noise, HG} (arb. units) ",
0790                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5);
0791   hMeanPedHGvsCellID->SetDirectory(0);
0792   TH2D* hspectraHGMeanVsLayer   = new TH2D( "hspectraHGMeanVsLayer","Mean Ped High Gain vs CellID; layer; brd channel; #mu_{Ped HG} (arb. units) ",
0793                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
0794   hspectraHGMeanVsLayer->SetDirectory(0);
0795   TH2D* hspectraHGSigmaVsLayer  = new TH2D( "hspectraHGSigmaVsLayer","Mean Ped High Gain vs CellID; layer; brd channel; #sigma_{Ped HG} (arb. units) ",
0796                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
0797   hspectraHGSigmaVsLayer->SetDirectory(0);
0798   TH1D* hMeanPedLGvsCellID      = new TH1D( "hMeanPedLG_vsCellID","mean Ped Low Gain vs CellID ; cell ID; #mu_{noise, LG} (arb. units) ",
0799                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5);
0800   hMeanPedLGvsCellID->SetDirectory(0);
0801   TH2D* hspectraLGMeanVsLayer   = new TH2D( "hspectraLGMeanVsLayer","Mean Ped Low Gain vs CellID; layer; brd channel; #mu_{PED LG} (arb. units) ",
0802                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
0803   hspectraLGMeanVsLayer->SetDirectory(0);
0804   TH2D* hspectraLGSigmaVsLayer  = new TH2D( "hspectraLGSigmaVsLayer","Mean Ped Low Gain vs CellID; layer; brd channel; #sigma_{Ped LG} (arb. units)",
0805                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
0806   hspectraLGSigmaVsLayer->SetDirectory(0);
0807
0808   std::map<int,TileSpectra> hSpectra;
0809   std::map<int, TileSpectra>::iterator ithSpectra;
0810
0811   std::cout << "Additional Output with histos being created: " << RootOutputNameHist.Data() << std::endl;
0812   if(Overwrite){
0813     std::cout << "recreating file with hists" << std::endl;
0814     RootOutputHist = new TFile(RootOutputNameHist.Data(),"RECREATE");
0815   } else{
0816     std::cout << "newly creating file with hists" << std::endl;
0817     RootOutputHist = new TFile(RootOutputNameHist.Data(),"CREATE");
0818   }
0819   // entering histoOutput file
0820   RootOutputHist->mkdir("IndividualCells");
0821   RootOutputHist->cd("IndividualCells");
0822
0823   RootOutput->cd();
0824   // Event loop to fill histograms & output tree
0825   std::cout << "N max layers: " << setup->GetNMaxLayer() << "\t columns: " <<  setup->GetNMaxColumn() << "\t row: " << setup->GetNMaxRow() << "\t module: " <<  setup->GetNMaxModule() << std::endl;
0826   if(TcalibIn) TcalibIn->GetEntry(0);
0827   int evts=TdataIn->GetEntries();
0828   int runNr = -1;
0829   for(int i=0; i<evts; i++){
0830     TdataIn->GetEntry(i);
0831     if (i == 0)runNr = event.GetRunNumber();
0832     if (i%5000 == 0&& i > 0 && debug > 0) std::cout << "Reading " <<  i << "/" << evts << " events"<< std::endl;
0833     for(int j=0; j<event.GetNTiles(); j++){
0834       Caen* aTile=(Caen*)event.GetTile(j);
0835       if (i == 0 && debug > 2){
0836         std::cout << ((TString)setup->DecodeCellID(aTile->GetCellID())).Data() << std::endl;
0837       }
0838       ithSpectra=hSpectra.find(aTile->GetCellID());
0839       if(ithSpectra!=hSpectra.end()){
0840         ithSpectra->second.Fill(aTile->GetADCLow(),aTile->GetADCHigh());
0841       } else {
0842         RootOutputHist->cd("IndividualCells");
0843         hSpectra[aTile->GetCellID()]=TileSpectra("1stExtraction",aTile->GetCellID(),calib.GetTileCalib(aTile->GetCellID()),debug);
0844         hSpectra[aTile->GetCellID()].Fill(aTile->GetADCLow(),aTile->GetADCHigh());
0845         RootOutput->cd();
0846       }
0847       hspectraHGvsCellID->Fill(aTile->GetCellID(), aTile->GetADCHigh());
0848       hspectraLGvsCellID->Fill(aTile->GetCellID(), aTile->GetADCLow());
0849     }
0850     RootOutput->cd();
0851     TdataOut->Fill();
0852   }
0853
0854   // fit pedestal
0855   double* parameters=new double[8];
0856   bool isGood;
0857
0858   for(ithSpectra=hSpectra.begin(); ithSpectra!=hSpectra.end(); ++ithSpectra){
0859     if ( debug > 2) std::cout << ((TString)setup->DecodeCellID(ithSpectra->second.GetCellID())).Data() << std::endl;
0860     isGood=ithSpectra->second.FitNoise(parameters, yearData, false);
0861     hMeanPedHGvsCellID->SetBinContent(hMeanPedHGvsCellID->FindBin(ithSpectra->second.GetCellID()), parameters[4]);
0862     hMeanPedHGvsCellID->SetBinError  (hMeanPedHGvsCellID->FindBin(ithSpectra->second.GetCellID()), parameters[6]);
0863     hMeanPedLGvsCellID->SetBinContent(hMeanPedLGvsCellID->FindBin(ithSpectra->second.GetCellID()), parameters[0]);
0864     hMeanPedLGvsCellID->SetBinError  (hMeanPedLGvsCellID->FindBin(ithSpectra->second.GetCellID()), parameters[2]);
0865
0866     int layer     = setup->GetLayer(ithSpectra->second.GetCellID());
0867     int chInLayer = setup->GetChannelInLayer(ithSpectra->second.GetCellID());
0868
0869     hspectraHGMeanVsLayer->SetBinContent(hspectraHGMeanVsLayer->FindBin(layer,chInLayer), parameters[4]);
0870     hspectraHGMeanVsLayer->SetBinError(hspectraHGMeanVsLayer->FindBin(layer,chInLayer), parameters[5]);
0871     hspectraHGSigmaVsLayer->SetBinContent(hspectraHGSigmaVsLayer->FindBin(layer,chInLayer), parameters[6]);
0872     hspectraHGSigmaVsLayer->SetBinError(hspectraHGSigmaVsLayer->FindBin(layer,chInLayer), parameters[7]);
0873     hspectraLGMeanVsLayer->SetBinContent(hspectraLGMeanVsLayer->FindBin(layer,chInLayer), parameters[0]);
0874     hspectraLGMeanVsLayer->SetBinError(hspectraLGMeanVsLayer->FindBin(layer,chInLayer), parameters[1]);
0875     hspectraLGSigmaVsLayer->SetBinContent(hspectraLGSigmaVsLayer->FindBin(layer,chInLayer), parameters[2]);
0876     hspectraLGSigmaVsLayer->SetBinError(hspectraLGSigmaVsLayer->FindBin(layer,chInLayer), parameters[3]);
0877   }
0878
0879   RootOutput->cd();
0880   // write output tree (copy of raw data)
0881   TdataOut->Write();
0882   // write setup tree (copy of raw data)
0883   TsetupIn->CloneTree()->Write();
0884
0885   if (IsCalibSaveToFile()){
0886     TString fileCalibPrint = RootOutputName;
0887     fileCalibPrint         = fileCalibPrint.ReplaceAll(".root","_calib.txt");
0888     calib.PrintCalibToFile(fileCalibPrint);
0889   }
0890
0891   TcalibOut->Fill();
0892   TcalibOut->Write();
0893   RootOutput->Write();
0894   RootOutput->Close();
0895
0896   // entering histoOutput file
0897   //RootOutputHist->cd();
0898   //RootOutputHist->mkdir("IndividualCells");
0899   RootOutputHist->cd("IndividualCells");
0900   for(ithSpectra=hSpectra.begin(); ithSpectra!=hSpectra.end(); ++ithSpectra){
0901     ithSpectra->second.Write(true);
0902   }
0903   RootOutputHist->cd();
0904   hspectraHGvsCellID->Write();
0905   hMeanPedHGvsCellID->Write();
0906   hspectraLGvsCellID->Write();
0907   hMeanPedLGvsCellID->Write();
0908   hspectraHGMeanVsLayer->Write();
0909   hspectraLGMeanVsLayer->Write();
0910   hspectraHGSigmaVsLayer->Write();
0911   hspectraLGSigmaVsLayer->Write();
0912
0913   // fill calib tree & write it
0914   // close open root files
0915   RootOutputHist->Write();
0916   RootOutputHist->Close();
0917
0918   RootInput->Close();
0919
0920
0921   // Get run info object
0922   std::map<int,RunInfo>::iterator it=ri.find(runNr);
0923
0924   // create directory for plot output
0925   TString outputDirPlots = GetPlotOutputDir();
0926   gSystem->Exec("mkdir -p "+outputDirPlots);
0927
0928   //**********************************************************************
0929   // Create canvases for channel overview plotting
0930   //**********************************************************************
0931   Double_t textSizeRel = 0.035;
0932   StyleSettingsBasics("pdf");
0933   SetPlotStyle();
0934
0935   TCanvas* canvas2DCorr = new TCanvas("canvasCorrPlots","",0,0,1450,1200);  // gives the page size
0936   DefaultCancasSettings( canvas2DCorr, 0.08, 0.13, 0.02, 0.07);
0937   canvas2DCorr->SetLogz();
0938
0939   PlotSimple2D( canvas2DCorr, hspectraHGvsCellID, 300, setup->GetMaxCellID()+1, textSizeRel, Form("%s/HG_Noise.pdf", outputDirPlots.Data()), it->second, 5 );
0940   PlotSimple2D( canvas2DCorr, hspectraLGvsCellID, 300, setup->GetMaxCellID()+1, textSizeRel, Form("%s/LG_Noise.pdf", outputDirPlots.Data()), it->second, 5);
0941
0942   canvas2DCorr->SetLogz(0);
0943   PlotSimple2D( canvas2DCorr, hspectraHGMeanVsLayer, -10000, -10000, textSizeRel, Form("%s/HG_NoiseMean.pdf", outputDirPlots.Data()), it->second, 5, kFALSE, "colz");
0944   PlotSimple2D( canvas2DCorr, hspectraHGSigmaVsLayer,-10000, -10000, textSizeRel, Form("%s/HG_NoiseSigma.pdf", outputDirPlots.Data()), it->second, 5,  kFALSE, "colz");
0945
0946   canvas2DCorr->SetLogz(0);
0947   PlotSimple2D( canvas2DCorr, hspectraLGMeanVsLayer, -10000, -10000, textSizeRel, Form("%s/LG_NoiseMean.pdf", outputDirPlots.Data()), it->second, 5, kFALSE, "colz");
0948   PlotSimple2D( canvas2DCorr, hspectraLGSigmaVsLayer, -10000, -10000, textSizeRel, Form("%s/LG_NoiseSigma.pdf", outputDirPlots.Data()), it->second, 5, kFALSE, "colz");
0949
0950
0951   //***********************************************************************************************************
0952   //********************************* 8 Panel overview plot  **************************************************
0953   //***********************************************************************************************************
0954   //*****************************************************************
0955     // Test beam geometry (beam coming from viewer)
0956     //===========================================================
0957     //||    8 (4)    ||    7 (5)   ||    6 (6)   ||    5 (7)   ||  row 0
0958     //===========================================================
0959     //||    1 (0)    ||    2 (1)   ||    3 (2)   ||    4 (3)   ||  row 1
0960     //===========================================================
0961     //    col 0     col 1       col 2     col  3
0962     // rebuild pad geom in similar way (numbering -1)
0963   //*****************************************************************
0964   TCanvas* canvas8Panel;
0965   TPad* pad8Panel[8];
0966   Double_t topRCornerX[8];
0967   Double_t topRCornerY[8];
0968   Int_t textSizePixel = 30;
0969   Double_t relSize8P[8];
0970   CreateCanvasAndPadsFor8PannelTBPlot(canvas8Panel, pad8Panel,  topRCornerX, topRCornerY, relSize8P, textSizePixel);
0971
0972   for (Int_t l = 0; l < setup->GetNMaxLayer()+1; l++){
0973     PlotNoiseWithFitsFullLayer (canvas8Panel,pad8Panel, topRCornerX, topRCornerY, relSize8P, textSizePixel,
0974                                 hSpectra, setup, true, 0, 275, 1.2, l, 0,
0975                                 Form("%s/Noise_HG_Layer%02d.pdf" ,outputDirPlots.Data(), l), it->second);
0976     PlotNoiseWithFitsFullLayer (canvas8Panel,pad8Panel, topRCornerX, topRCornerY, relSize8P, textSizePixel,
0977                                 hSpectra, setup, false, 0, 275, 1.2, l, 0,
0978                                 Form("%s/Noise_LG_Layer%02d.pdf" ,outputDirPlots.Data(), l), it->second);
0979   }
0980
0981
0982   return true;
0983 }
0984
0985 // ****************************************************************************
0986 // extract pedestral from dedicated data run, no other data in run
0987 // ****************************************************************************
0988 bool Analyses::CorrectPedestal(void){
0989   std::cout<<"Correct Pedestal"<<std::endl;
0990   int evts=TdataIn->GetEntries();
0991
0992   std::map<int,RunInfo> ri=readRunInfosFromFile(RunListInputName.Data(),debug,0);
0993
0994   std::map<int,TileSpectra> hSpectra;
0995   std::map<int, TileSpectra>::iterator ithSpectra;
0996   TcalibIn->GetEntry(0);
0997   int runNr = -1;
0998
0999   std::map<int,short> bcmap;
1000   std::map<int,short>::iterator bcmapIte;
1001   if (CalcBadChannel == 1)
1002     bcmap = ReadExternalBadChannelMap();
1003
1004   ROOT::Math::MinimizerOptions::SetDefaultMinimizer("Minuit2", "Migrad");
1005
1006   for(int i=0; i<evts; i++){
1007     if (i%5000 == 0&& i > 0 && debug > 0) std::cout << "Reading " <<  i << "/" << evts << " events"<< std::endl;
1008     TdataIn->GetEntry(i);
1009     if (i == 0)runNr = event.GetRunNumber();
1010
1011     if (CalcBadChannel > 0 || ExtPlot > 0){
1012       for(int j=0; j<event.GetNTiles(); j++){
1013         Caen* aTile=(Caen*)event.GetTile(j);
1014         ithSpectra=hSpectra.find(aTile->GetCellID());
1015         double hgCorr = aTile->GetADCHigh()-calib.GetPedestalMeanH(aTile->GetCellID());
1016         double lgCorr = aTile->GetADCLow()-calib.GetPedestalMeanL(aTile->GetCellID());
1017
1018         if(ithSpectra!=hSpectra.end()){
1019           ithSpectra->second.FillSpectra(lgCorr,hgCorr);
1020         } else {
1021           hSpectra[aTile->GetCellID()]=TileSpectra("ped",aTile->GetCellID(),calib.GetTileCalib(aTile->GetCellID()),debug);
1022           hSpectra[aTile->GetCellID()].FillSpectra(lgCorr,hgCorr);;
1023         }
1024       }
1025     }
1026     RootOutput->cd();
1027     TdataOut->Fill();
1028   }
1029   //RootOutput->cd();
1030   TdataOut->Write();
1031   TsetupIn->CloneTree()->Write();
1032
1033   std::map<int,RunInfo>::iterator it=ri.find(runNr);
1034   TString outputDirPlots = GetPlotOutputDir();
1035   Double_t textSizeRel = 0.035;
1036
1037   if (CalcBadChannel > 0 || ExtPlot > 0) {
1038     gSystem->Exec("mkdir -p "+outputDirPlots);
1039     StyleSettingsBasics("pdf");
1040     SetPlotStyle();
1041   }
1042
1043   if (CalcBadChannel > 0 ){
1044     //***********************************************************************************************
1045     //***** Monitoring histos for calib parameters & fits results of 1st iteration ******************
1046     //***********************************************************************************************
1047     int maxChannelPerLayer        = (setup->GetNMaxColumn()+1)*(setup->GetNMaxRow()+1);
1048     // monitoring applied pedestals
1049     TH1D* hBCvsCellID      = new TH1D( "hBC_vsCellID","Bad Channel vs CellID ; cell ID; BC flag ",
1050                                               setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5);
1051     hBCvsCellID->SetDirectory(0);
1052     TH2D* hBCVsLayer   = new TH2D( "hBCVsLayer","Bad Channel Map; layer; brd channel; BC flag  ",
1053                                               setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
1054     hBCVsLayer->SetDirectory(0);
1055
1056     int currCells = 0;
1057     if ( debug > 0 && CalcBadChannel == 2)
1058       std::cout << "============================== starting bad channel extraction" << std::endl;
1059     if ( debug > 0 && CalcBadChannel == 1)
1060       std::cout << "============================== setting bad channel according to external map" << std::endl;
1061
1062     for(ithSpectra=hSpectra.begin(); ithSpectra!=hSpectra.end(); ++ithSpectra){
1063       if (currCells%20 == 0 && currCells > 0 && debug > 0)
1064         std::cout << "============================== cell " <<  currCells << " / " << hSpectra.size() << " cells" << std::endl;
1065       currCells++;
1066       short bc   = 3;
1067       if (CalcBadChannel == 2){
1068         bc        = ithSpectra->second.DetermineBadChannel();
1069       } else if (CalcBadChannel == 1 && bcmap.size() > 0 ) {
1070         bcmapIte  = bcmap.find(ithSpectra->second.GetCellID());
1071         if(bcmapIte!=bcmap.end())
1072           bc        = bcmapIte->second;
1073         else
1074           bc        = 3;
1075       }
1076
1077       long cellID   = ithSpectra->second.GetCellID();
1078       if (CalcBadChannel == 1)
1079         ithSpectra->second.SetBadChannelInCalib(bc);
1080
1081       int layer     = setup->GetLayer(cellID);
1082       int chInLayer = setup->GetChannelInLayer(cellID);
1083       if (debug > 0 && bc > -1 && bc < 3)
1084         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;
1085
1086       hBCvsCellID->SetBinContent(hBCvsCellID->FindBin(cellID), calib.GetBadChannel(cellID));
1087       int bin2D     = hBCVsLayer->FindBin(layer,chInLayer);
1088       hBCVsLayer->SetBinContent(bin2D, calib.GetBadChannel(cellID));
1089     }
1090     hBCvsCellID->Write();
1091     hBCVsLayer->Write();
1092
1093     //**********************************************************************
1094     // Create canvases for channel overview plotting
1095     //**********************************************************************
1096     // Get run info object
1097     // create directory for plot output
1098
1099     TCanvas* canvas2DCorr = new TCanvas("canvasCorrPlots","",0,0,1450,1300);  // gives the page size
1100     DefaultCancasSettings( canvas2DCorr, 0.08, 0.13, 0.045, 0.07);
1101
1102     canvas2DCorr->SetLogz(0);
1103
1104     PlotSimple2DZRange( canvas2DCorr, hBCVsLayer, -10000, -10000, -0.1, 3.1, textSizeRel, Form("%s/BadChannelMap.pdf", outputDirPlots.Data()), it->second, 1, "colz", true);
1105     calib.SetBCCalib(true);
1106   }
1107
1108   if (ExtPlot > 0){
1109     //***********************************************************************************************************
1110     //********************************* 8 Panel overview plot  **************************************************
1111     //***********************************************************************************************************
1112     //*****************************************************************
1113       // Test beam geometry (beam coming from viewer)
1114       //===========================================================
1115       //||    8 (4)    ||    7 (5)   ||    6 (6)   ||    5 (7)   ||  row 0
1116       //===========================================================
1117       //||    1 (0)    ||    2 (1)   ||    3 (2)   ||    4 (3)   ||  row 1
1118       //===========================================================
1119       //    col 0     col 1       col 2     col  3
1120       // rebuild pad geom in similar way (numbering -1)
1121     //*****************************************************************
1122     TCanvas* canvas8Panel;
1123     TPad* pad8Panel[8];
1124     Double_t topRCornerX[8];
1125     Double_t topRCornerY[8];
1126     Int_t textSizePixel = 30;
1127     Double_t relSize8P[8];
1128     CreateCanvasAndPadsFor8PannelTBPlot(canvas8Panel, pad8Panel,  topRCornerX, topRCornerY, relSize8P, textSizePixel);
1129
1130     calib.PrintGlobalInfo();
1131     // Double_t maxHG = ReturnMipPlotRangeDepVov(calib.GetVov(),true);
1132     // Double_t maxLG = ReturnMipPlotRangeDepVov(calib.GetVov(),false);
1133     Double_t maxHG = 600;
1134     Double_t maxLG = 200;
1135
1136     std::cout << "plotting single layers" << std::endl;
1137     for (Int_t l = 0; l < setup->GetNMaxLayer()+1; l++){
1138       if (l%10 == 0 && l > 0 && debug > 0)
1139         std::cout << "============================== layer " <<  l << " / " << setup->GetNMaxLayer() << " layers" << std::endl;
1140       PlotNoiseWithFitsFullLayer (canvas8Panel,pad8Panel, topRCornerX, topRCornerY, relSize8P, textSizePixel,
1141                                   hSpectra, setup, true, -100, maxHG, 1.2, l, 0,
1142                                   Form("%s/SpectraWithNoiseFit_HG_Layer%02d.pdf" ,outputDirPlots.Data(), l), it->second);
1143       PlotNoiseWithFitsFullLayer (canvas8Panel,pad8Panel, topRCornerX, topRCornerY, relSize8P, textSizePixel,
1144                                   hSpectra, setup, false, -20, maxLG, 1.2, l, 0,
1145                                   Form("%s/SpectraWithNoiseFit_LG_Layer%02d.pdf" ,outputDirPlots.Data(), l), it->second);
1146     }
1147     std::cout << "ending plotting single layers" << std::endl;
1148   }
1149
1150   RootOutput->cd();
1151
1152   std::cout<<"What is the value? <ped mean high>: "<<calib.GetAveragePedestalMeanHigh() << "\t good channels: " << calib.GetNumberOfChannelsWithBCflag(3) <<std::endl;
1153   if (IsCalibSaveToFile()){
1154     TString fileCalibPrint = RootOutputName;
1155     fileCalibPrint         = fileCalibPrint.ReplaceAll(".root","_calib.txt");
1156     calib.PrintCalibToFile(fileCalibPrint);
1157   }
1158
1159   TcalibOut->Fill();
1160   TcalibOut->Write();
1161
1162   RootOutput->Close();
1163   RootInput->Close();
1164   return true;
1165 }
1166
1167
1168 //***********************************************************************************************
1169 //*********************** intial scaling calculation function *********************************
1170 //***********************************************************************************************
1171 bool Analyses::GetScaling(void){
1172   std::cout<<"GetScaling"<<std::endl;
1173
1174   std::map<int,RunInfo> ri=readRunInfosFromFile(RunListInputName.Data(),debug,0);
1175
1176   std::map<int,TileSpectra> hSpectra;
1177   std::map<int,TileSpectra> hSpectraTrigg;
1178   std::map<int, TileSpectra>::iterator ithSpectra;
1179   std::map<int, TileSpectra>::iterator ithSpectraTrigg;
1180
1181   ROOT::Math::MinimizerOptions::SetDefaultMinimizer("Minuit2", "Migrad");
1182
1183   std::cout << "Additional Output with histos being created: " << RootOutputNameHist.Data() << std::endl;
1184   if(Overwrite){
1185     std::cout << "recreating file with hists" << std::endl;
1186     RootOutputHist = new TFile(RootOutputNameHist.Data(),"RECREATE");
1187   } else{
1188     std::cout << "newly creating file with hists" << std::endl;
1189     RootOutputHist = new TFile(RootOutputNameHist.Data(),"CREATE");
1190   }
1191
1192   //***********************************************************************************************
1193   //************************* first pass over tree to extract spectra *****************************
1194   //***********************************************************************************************
1195   // entering histoOutput file
1196   RootOutputHist->mkdir("IndividualCells");
1197   RootOutputHist->cd("IndividualCells");
1198
1199   TcalibIn->GetEntry(0);
1200   int evts=TdataIn->GetEntries();
1201   int runNr = -1;
1202   for(int i=0; i<evts; i++){
1203     TdataIn->GetEntry(i);
1204     if (i == 0)runNr = event.GetRunNumber();
1205     if (i%5000 == 0 && i > 0 && debug > 0) std::cout << "Reading " <<  i << " / " << evts << " events" << std::endl;
1206     if (i == 0 && debug > 2) std::cout << "decoding cell IDs" << std::endl ;
1207     for(int j=0; j<event.GetNTiles(); j++){
1208       Caen* aTile=(Caen*)event.GetTile(j);
1209       if (i == 0 && debug > 2) std::cout << ((TString)setup->DecodeCellID(aTile->GetCellID())).Data() << std::endl;
1210
1211       ithSpectra=hSpectra.find(aTile->GetCellID());
1212       double hgCorr = aTile->GetADCHigh()-calib.GetPedestalMeanH(aTile->GetCellID());
1213       double lgCorr = aTile->GetADCLow()-calib.GetPedestalMeanL(aTile->GetCellID());
1214
1215       if(ithSpectra!=hSpectra.end()){
1216         ithSpectra->second.FillSpectra(lgCorr,hgCorr);
1217         if (hgCorr > 3*calib.GetPedestalSigH(aTile->GetCellID()) && lgCorr > 3*calib.GetPedestalSigL(aTile->GetCellID()) && hgCorr < 3900 )
1218           ithSpectra->second.FillCorr(lgCorr,hgCorr);
1219       } else {
1220         RootOutputHist->cd("IndividualCells");
1221         hSpectra[aTile->GetCellID()]=TileSpectra("mip1st",aTile->GetCellID(),calib.GetTileCalib(aTile->GetCellID()),debug);
1222         hSpectra[aTile->GetCellID()].FillSpectra(lgCorr,hgCorr);;
1223         if (hgCorr > 3*calib.GetPedestalSigH(aTile->GetCellID()) && lgCorr > 3*calib.GetPedestalSigL(aTile->GetCellID() && hgCorr < 3900) )
1224           hSpectra[aTile->GetCellID()].FillCorr(lgCorr,hgCorr);;
1225         RootOutput->cd();
1226       }
1227     }
1228     RootOutput->cd();
1229   }
1230   // TdataOut->Write();
1231   TsetupIn->CloneTree()->Write();
1232
1233   RootOutputHist->cd();
1234
1235   //***********************************************************************************************
1236   //***** Monitoring histos for calib parameters & fits results of 1st iteration ******************
1237   //***********************************************************************************************
1238   int maxChannelPerLayer        = (setup->GetNMaxColumn()+1)*(setup->GetNMaxRow()+1);
1239   // monitoring applied pedestals
1240   TH1D* hMeanPedHGvsCellID      = new TH1D( "hMeanPedHG_vsCellID","mean Ped High Gain vs CellID ; cell ID; #mu_{noise, HG} (arb. units) ",
1241                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5);
1242   hMeanPedHGvsCellID->SetDirectory(0);
1243   TH1D* hMeanPedHG              = new TH1D( "hMeanPedHG","mean Ped High Gain ; #mu_{noise, HG} (arb. units); counts ",
1244                                             500, -0.5, 500-0.5);
1245   hMeanPedHG->SetDirectory(0);
1246   TH2D* hspectraHGMeanVsLayer   = new TH2D( "hspectraHGMeanVsLayer","Mean Ped High Gain; layer; brd channel; #mu_{Ped HG} (arb. units) ",
1247                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
1248   hspectraHGMeanVsLayer->SetDirectory(0);
1249   TH2D* hspectraHGSigmaVsLayer  = new TH2D( "hspectraHGSigmaVsLayer","Sigma Ped High Gain; layer; brd channel; #sigma_{Ped HG} (arb. units) ",
1250                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
1251   hspectraHGSigmaVsLayer->SetDirectory(0);
1252   TH1D* hMeanPedLGvsCellID      = new TH1D( "hMeanPedLG_vsCellID","mean Ped Low Gain vs CellID ; cell ID; #mu_{noise, LG} (arb. units) ",
1253                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5);
1254   hMeanPedLGvsCellID->SetDirectory(0);
1255   TH1D* hMeanPedLG             = new TH1D( "hMeanPedLG","mean Ped Low Gain ; #mu_{noise, LG} (arb. units); counts ",
1256                                             500, -0.5, 500-0.5);
1257   hMeanPedLG->SetDirectory(0);
1258   TH2D* hspectraLGMeanVsLayer   = new TH2D( "hspectraLGMeanVsLayer","Mean Ped Low Gain; layer; brd channel; #mu_{PED LG} (arb. units) ",
1259                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
1260   hspectraLGMeanVsLayer->SetDirectory(0);
1261   TH2D* hspectraLGSigmaVsLayer  = new TH2D( "hspectraLGSigmaVsLayer","Sigma Ped Low Gain; layer; brd channel; #sigma_{Ped LG} (arb. units)",
1262                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
1263   hspectraLGSigmaVsLayer->SetDirectory(0);
1264   // monitoring 1st iteration mip fits
1265   TH2D* hspectraHGMaxVsLayer1st   = new TH2D( "hspectraHGMaxVsLayer_1st","Max High Gain; layer; brd channel; Max_{HG, 1^{st}} (arb. units) ",
1266                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
1267   hspectraHGMaxVsLayer1st->SetDirectory(0);
1268   TH2D* hspectraHGFWHMVsLayer1st   = new TH2D( "hspectraHGFWHMVsLayer_1st","FWHM High Gain; layer; brd channel; FWHM_{HG, 1^{st}} (arb. units) ",
1269                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
1270   hspectraHGFWHMVsLayer1st->SetDirectory(0);
1271   TH2D* hspectraHGLMPVVsLayer1st   = new TH2D( "hspectraHGLMPVVsLayer_1st","MPV High Gain; layer; brd channel; MPV_{HG, 1^{st}} (arb. units) ",
1272                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
1273   hspectraHGLMPVVsLayer1st->SetDirectory(0);
1274   TH2D* hspectraHGLSigmaVsLayer1st = new TH2D( "hspectraHGLSigmaVsLayer_1st","Sigma Landau High Gain; layer; brd channel; #sigma_{L,HG, 1^{st}} (arb. units) ",
1275                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
1276   hspectraHGLSigmaVsLayer1st->SetDirectory(0);
1277   TH2D* hspectraHGGSigmaVsLayer1st = new TH2D( "hspectraHGGSigmaVsLayer_1st","Sigma Gauss High Gain; layer; brd channel; #sigma_{G,HG, 1^{st}} (arb. units) ",
1278                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
1279   hspectraHGGSigmaVsLayer1st->SetDirectory(0);
1280   TH2D* hLGHGCorrVsLayer = new TH2D( "hLGHGCorrVsLayer","LG-HG corr; layer; brd channel; a_{LG-HG} (arb. units) ",
1281                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
1282   hLGHGCorrVsLayer->SetDirectory(0);
1283   TH2D* hHGLGCorrVsLayer = new TH2D( "hHGLGCorrVsLayer","HG-LG corr; layer; brd channel; a_{HG-LG} (arb. units) ",
1284                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
1285   hHGLGCorrVsLayer->SetDirectory(0);
1286
1287   TH1D* hMaxHG1st             = new TH1D( "hMaxHG1st","Max High Gain ;Max_{HG, 1^{st}} (arb. units) ; counts ",
1288                                             2000, -0.5, 2000-0.5);
1289   hMaxHG1st->SetDirectory(0);
1290   TH1D* hLGHGCorr             = new TH1D( "hLGHGCorr","LG-HG corr ; a_{LG-HG} (arb. units) ; counts ",
1291                                             400, -20, 20);
1292   hLGHGCorr->SetDirectory(0);
1293   TH1D* hHGLGCorr             = new TH1D( "hHGLGCorr","LG-HG corr ; a_{HG-LG} (arb. units) ; counts ",
1294                                             400, -1., 1.);
1295   hHGLGCorr->SetDirectory(0);
1296
1297
1298   // fit pedestal
1299   double* parameters    = new double[6];
1300   double* parErrAndRes  = new double[6];
1301   bool isGood;
1302   int currCells = 0;
1303   if ( debug > 0)
1304     std::cout << "============================== starting fitting 1st iteration" << std::endl;
1305
1306   for(ithSpectra=hSpectra.begin(); ithSpectra!=hSpectra.end(); ++ithSpectra){
1307     if (currCells%20 == 0 && currCells > 0 && debug > 0)
1308       std::cout << "============================== cell " <<  currCells << " / " << hSpectra.size() << " cells" << std::endl;
1309     currCells++;
1310     for (int p = 0; p < 6; p++){
1311       parameters[p]   = 0;
1312       parErrAndRes[p] = 0;
1313     }
1314     isGood        = ithSpectra->second.FitMipHG(parameters, parErrAndRes, debug, yearData, false, calib.GetVov(), 1);
1315     long cellID   = ithSpectra->second.GetCellID();
1316     int layer     = setup->GetLayer(cellID);
1317     int chInLayer = setup->GetChannelInLayer(cellID);
1318
1319     // fill cross-check histos
1320     hMeanPedHGvsCellID->SetBinContent(hMeanPedHGvsCellID->FindBin(cellID), calib.GetPedestalMeanH(cellID));
1321     hMeanPedLGvsCellID->SetBinContent(hMeanPedLGvsCellID->FindBin(cellID), calib.GetPedestalMeanL(cellID));
1322     hMeanPedHG->Fill(calib.GetPedestalMeanH(cellID));
1323     hMeanPedLG->Fill(calib.GetPedestalMeanL(cellID));
1324
1325     int bin2D     = hspectraHGMeanVsLayer->FindBin(layer,chInLayer);
1326     hspectraHGMeanVsLayer->SetBinContent(bin2D, calib.GetPedestalMeanH(cellID));
1327     hspectraHGSigmaVsLayer->SetBinContent(bin2D, calib.GetPedestalSigH(cellID));
1328     hspectraLGMeanVsLayer->SetBinContent(bin2D, calib.GetPedestalMeanL(cellID));
1329     hspectraLGSigmaVsLayer->SetBinContent(bin2D, calib.GetPedestalSigL(cellID));
1330
1331     if (isGood){
1332       hspectraHGMaxVsLayer1st->SetBinContent(bin2D, parameters[4]);
1333       hspectraHGFWHMVsLayer1st->SetBinContent(bin2D, parameters[5]);
1334       hspectraHGLMPVVsLayer1st->SetBinContent(bin2D, parameters[1]);
1335       hspectraHGLMPVVsLayer1st->SetBinError(bin2D, parErrAndRes[1]);
1336       hspectraHGLSigmaVsLayer1st->SetBinContent(bin2D, parameters[0]);
1337       hspectraHGLSigmaVsLayer1st->SetBinError(bin2D, parErrAndRes[0]);
1338       hspectraHGGSigmaVsLayer1st->SetBinContent(bin2D, parameters[3]);
1339       hspectraHGGSigmaVsLayer1st->SetBinError(bin2D, parErrAndRes[3]);
1340
1341       hMaxHG1st->Fill(parameters[4]);
1342     }
1343
1344     isGood=ithSpectra->second.FitCorr(debug);
1345     if (ithSpectra->second.GetCorrModel(0)){
1346       hLGHGCorrVsLayer->SetBinContent(bin2D,ithSpectra->second.GetCorrModel(0)->GetParameter(1));
1347       hLGHGCorrVsLayer->SetBinError(bin2D,ithSpectra->second.GetCorrModel(0)->GetParError(1));
1348       hLGHGCorr->Fill(ithSpectra->second.GetCorrModel(0)->GetParameter(1));
1349     }
1350     if (ithSpectra->second.GetCorrModel(1)){
1351       hHGLGCorrVsLayer->SetBinContent(bin2D,ithSpectra->second.GetCorrModel(1)->GetParameter(1));
1352       hHGLGCorrVsLayer->SetBinError(bin2D,ithSpectra->second.GetCorrModel(1)->GetParError(1));
1353       hHGLGCorr->Fill(ithSpectra->second.GetCorrModel(1)->GetParameter(1));
1354     }
1355   }
1356   if ( debug > 0)
1357     std::cout << "============================== done fitting 1st iteration" << std::endl;
1358
1359   TH1D* hHGTileSum[20];
1360   for (int c = 0; c < maxChannelPerLayer; c++ ){
1361     hHGTileSum[c] = new TH1D( Form("hHGTileSum_absChB%d",c),"av ADC surrounding cells ; ADC (arb. units); counts ",
1362                               4000, -0.5, 4000-0.5);
1363     hHGTileSum[c]->SetDirectory(0);
1364   }
1365
1366   // setup trigger sel
1367   TRandom3* rand    = new TRandom3();
1368   Int_t localTriggerTiles = 4;
1369   double factorMinTrigg   = 0.5;
1370   double factorMaxTrigg   = 4.;
1371   if (yearData == 2023){
1372     localTriggerTiles = 6;
1373     factorMaxTrigg    = 3.;
1374   }
1375   RootOutputHist->mkdir("IndividualCellsTrigg");
1376   RootOutputHist->cd("IndividualCellsTrigg");
1377   //***********************************************************************************************
1378   //************************* first pass over tree to extract spectra *****************************
1379   //***********************************************************************************************
1380   int actCh1st = 0;
1381   double averageScale = calib.GetAverageScaleHigh(actCh1st);
1382   double meanFWHMHG   = calib.GetAverageScaleWidthHigh();
1383   double avHGLGCorr   = calib.GetAverageHGLGCorr();
1384   double avLGHGCorr   = calib.GetAverageLGHGCorr();
1385   std::cout << "average HG mip: " << averageScale << "\t active ch: "<< actCh1st<< std::endl;
1386   for(int i=0; i<evts; i++){
1387     TdataIn->GetEntry(i);
1388     if (i%5000 == 0 && i > 0 && debug > 0) std::cout << "Reading " <<  i << " / " << evts << " events" << std::endl;
1389     for(int j=0; j<event.GetNTiles(); j++){
1390       Caen* aTile=(Caen*)event.GetTile(j);
1391       if (i == 0 && debug > 2) std::cout << ((TString)setup->DecodeCellID(aTile->GetCellID())).Data() << std::endl;
1392       long currCellID = aTile->GetCellID();
1393
1394       // read current tile
1395       ithSpectraTrigg=hSpectraTrigg.find(aTile->GetCellID());
1396       double hgCorr = aTile->GetADCHigh()-calib.GetPedestalMeanH(aTile->GetCellID());
1397       double lgCorr = aTile->GetADCLow()-calib.GetPedestalMeanL(aTile->GetCellID());
1398
1399       aTile->SetLocalTriggerPrimitive(event.CalculateLocalMuonTrigg(calib, rand, currCellID, localTriggerTiles, avLGHGCorr));
1400       // estimate local muon trigger
1401       bool localMuonTrigg = event.InspectIfLocalMuonTrigg(currCellID, averageScale, factorMinTrigg, factorMaxTrigg);
1402       int chInLayer = setup->GetChannelInLayer(currCellID);
1403       hHGTileSum[chInLayer]->Fill(aTile->GetLocalTriggerPrimitive());
1404
1405       if(ithSpectraTrigg!=hSpectraTrigg.end()){
1406         ithSpectraTrigg->second.FillTrigger(aTile->GetLocalTriggerPrimitive());
1407       } else {
1408         RootOutputHist->cd("IndividualCellsTrigg");
1409         hSpectraTrigg[currCellID]=TileSpectra("mipTrigg",currCellID,calib.GetTileCalib(currCellID),debug);
1410         hSpectraTrigg[currCellID].FillTrigger(aTile->GetLocalTriggerPrimitive());;
1411         RootOutput->cd();
1412       }
1413
1414       // only fill tile spectra if 4 surrounding tiles on average are compatible with muon response
1415       if (localMuonTrigg){
1416         aTile->SetLocalTriggerBit(1);
1417         ithSpectraTrigg=hSpectraTrigg.find(aTile->GetCellID());
1418         ithSpectraTrigg->second.FillSpectra(lgCorr,hgCorr);
1419         if (hgCorr > 3*calib.GetPedestalSigH(currCellID) && lgCorr > 3*calib.GetPedestalSigL(currCellID) && hgCorr < 3900 )
1420           ithSpectraTrigg->second.FillCorr(lgCorr,hgCorr);
1421       }
1422     }
1423     RootOutput->cd();
1424     TdataOut->Fill();
1425   }
1426   TdataOut->Write();
1427
1428   //***********************************************************************************************
1429   //***** Monitoring histos for fits results of 2nd iteration ******************
1430   //***********************************************************************************************
1431   RootOutputHist->cd();
1432
1433   // monitoring trigger
1434   TH2D* hmipTriggers              = new TH2D( "hmipTriggers","muon triggers; layer; brd channel; counts ",
1435                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
1436   hmipTriggers->SetDirectory(0);
1437   TH2D* hSuppresionNoise          = new TH2D( "hSuppresionNoise","S/B noise region; layer; brd channel; S/B noise region",
1438                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
1439   hSuppresionNoise->SetDirectory(0);
1440   TH2D* hSuppresionSignal         = new TH2D( "hSuppresionSignal","S/B signal region; layer; brd channel; S/B signal region",
1441                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
1442   hSuppresionSignal->SetDirectory(0);
1443
1444   // monitoring 2nd iteration mip fits
1445   TH2D* hspectraHGMaxVsLayer2nd   = new TH2D( "hspectraHGMaxVsLayer_2nd","Max High Gain; layer; brd channel; Max_{HG, 2^{nd}} (arb. units) ",
1446                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
1447   hspectraHGMaxVsLayer2nd->SetDirectory(0);
1448   TH2D* hspectraHGFWHMVsLayer2nd   = new TH2D( "hspectraHGFWHMVsLayer_2nd","FWHM High Gain; layer; brd channel; FWHM_{HG, 2^{nd}} (arb. units) ",
1449                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
1450   hspectraHGFWHMVsLayer2nd->SetDirectory(0);
1451   TH2D* hspectraHGLMPVVsLayer2nd   = new TH2D( "hspectraHGLMPVVsLayer_2nd","MPV High Gain; layer; brd channel; MPV_{HG, 2^{nd}} (arb. units) ",
1452                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
1453   hspectraHGLMPVVsLayer2nd->SetDirectory(0);
1454   TH2D* hspectraHGLSigmaVsLayer2nd = new TH2D( "hspectraHGLSigmaVsLayer_2nd","Sigma Landau High Gain; layer; brd channel; #sigma_{L,HG, 2^{nd}} (arb. units) ",
1455                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
1456   hspectraHGLSigmaVsLayer2nd->SetDirectory(0);
1457   TH2D* hspectraHGGSigmaVsLayer2nd = new TH2D( "hspectraHGGSigmaVsLayer_2nd","Sigma Gauss High Gain; layer; brd channel; #sigma_{G,HG, 2^{nd}} (arb. units) ",
1458                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
1459   hspectraHGGSigmaVsLayer2nd->SetDirectory(0);
1460   TH2D* hspectraLGMaxVsLayer2nd   = new TH2D( "hspectraLGMaxVsLayer_2nd","Max High Gain; layer; brd channel; Max_{LG, 2^{nd}} (arb. units) ",
1461                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
1462   hspectraLGMaxVsLayer2nd->SetDirectory(0);
1463   TH2D* hspectraLGFWHMVsLayer2nd   = new TH2D( "hspectraLGFWHMVsLayer_2nd","FWHM High Gain; layer; brd channel; FWHM_{LG, 2^{nd}} (arb. units) ",
1464                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
1465   hspectraLGFWHMVsLayer2nd->SetDirectory(0);
1466   TH2D* hspectraLGLMPVVsLayer2nd   = new TH2D( "hspectraLGLMPVVsLayer_2nd","MPV High Gain; layer; brd channel; MPV_{LG, 2^{nd}} (arb. units) ",
1467                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
1468   hspectraLGLMPVVsLayer2nd->SetDirectory(0);
1469   TH2D* hspectraLGLSigmaVsLayer2nd = new TH2D( "hspectraLGLSigmaVsLayer_2nd","Sigma Landau High Gain; layer; brd channel; #sigma_{L,LG, 2^{nd}} (arb. units) ",
1470                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
1471   hspectraLGLSigmaVsLayer2nd->SetDirectory(0);
1472   TH2D* hspectraLGGSigmaVsLayer2nd = new TH2D( "hspectraLGGSigmaVsLayer_2nd","Sigma Gauss High Gain; layer; brd channel; #sigma_{G,LG, 2^{nd}} (arb. units) ",
1473                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
1474   hspectraLGGSigmaVsLayer2nd->SetDirectory(0);
1475
1476   TH1D* hMaxHG2nd             = new TH1D( "hMaxHG2nd","Max High Gain ;Max_{HG, 2^{nd}} (arb. units) ; counts ",
1477                                             2000, -0.5, 2000-0.5);
1478   hMaxHG2nd->SetDirectory(0);
1479   TH1D* hMaxLG2nd             = new TH1D( "hMaxLG2nd","Max High Gain ;Max_{LG, 2^{nd}} (arb. units) ; counts ",
1480                                             400, -0.5, 400-0.5);
1481   hMaxLG2nd->SetDirectory(0);
1482
1483
1484   currCells = 0;
1485   if ( debug > 0)
1486     std::cout << "============================== starting fitting 2nd iteration" << std::endl;
1487   for(ithSpectraTrigg=hSpectraTrigg.begin(); ithSpectraTrigg!=hSpectraTrigg.end(); ++ithSpectraTrigg){
1488     if (currCells%20 == 0 && currCells > 0 && debug > 0)
1489       std::cout << "============================== cell " <<  currCells << " / " << hSpectraTrigg.size() << " cells" << std::endl;
1490     currCells++;
1491     for (int p = 0; p < 6; p++){
1492       parameters[p]   = 0;
1493       parErrAndRes[p] = 0;
1494     }
1495     isGood=ithSpectraTrigg->second.FitMipHG(parameters, parErrAndRes, debug, yearData, false, calib.GetVov(), averageScale);
1496
1497     long cellID     = ithSpectraTrigg->second.GetCellID();
1498     int layer       = setup->GetLayer(cellID);
1499     int chInLayer   = setup->GetChannelInLayer(cellID);
1500     int bin2D       = hspectraHGMeanVsLayer->FindBin(layer,chInLayer);
1501
1502     Int_t binNLow   = ithSpectraTrigg->second.GetHG()->FindBin(-1*calib.GetPedestalSigH(cellID));
1503     Int_t binNHigh  = ithSpectraTrigg->second.GetHG()->FindBin(3*calib.GetPedestalSigH(cellID));
1504     Int_t binSHigh  = ithSpectraTrigg->second.GetHG()->FindBin(3800);
1505
1506     double S_NoiseR = ithSpectraTrigg->second.GetHG()->Integral(binNLow, binNHigh);
1507     double S_SigR   = ithSpectraTrigg->second.GetHG()->Integral(binNHigh, binSHigh);
1508
1509     ithSpectra      = hSpectra.find(cellID);
1510     double B_NoiseR = ithSpectra->second.GetHG()->Integral(binNLow , binNHigh);
1511     double B_SigR   = ithSpectra->second.GetHG()->Integral(binNHigh, binSHigh);
1512
1513     double SB_NoiseR  = (B_NoiseR != 0.) ? S_NoiseR/B_NoiseR : 0;
1514     double SB_SigR    = (B_SigR != 0.) ? S_SigR/B_SigR : 0;
1515
1516     hmipTriggers->SetBinContent(bin2D, ithSpectraTrigg->second.GetHG()->GetEntries());
1517     hSuppresionNoise->SetBinContent(bin2D, SB_NoiseR);
1518     hSuppresionSignal->SetBinContent(bin2D, SB_SigR);
1519     if (isGood){
1520       hspectraHGMaxVsLayer2nd->SetBinContent(bin2D, parameters[4]);
1521       hspectraHGFWHMVsLayer2nd->SetBinContent(bin2D, parameters[5]);
1522       hspectraHGLMPVVsLayer2nd->SetBinContent(bin2D, parameters[1]);
1523       hspectraHGLMPVVsLayer2nd->SetBinError(bin2D, parErrAndRes[1]);
1524       hspectraHGLSigmaVsLayer2nd->SetBinContent(bin2D, parameters[0]);
1525       hspectraHGLSigmaVsLayer2nd->SetBinError(bin2D, parErrAndRes[0]);
1526       hspectraHGGSigmaVsLayer2nd->SetBinContent(bin2D, parameters[3]);
1527       hspectraHGGSigmaVsLayer2nd->SetBinError(bin2D, parErrAndRes[3]);
1528       hMaxHG2nd->Fill(parameters[4]);
1529     }
1530     for (int p = 0; p < 6; p++){
1531       parameters[p]   = 0;
1532       parErrAndRes[p] = 0;
1533     }
1534     isGood=ithSpectraTrigg->second.FitMipLG(parameters, parErrAndRes, debug, yearData, false, 1);
1535     if (isGood){
1536       hspectraLGMaxVsLayer2nd->SetBinContent(bin2D, parameters[4]);
1537       hspectraLGFWHMVsLayer2nd->SetBinContent(bin2D, parameters[5]);
1538       hspectraLGLMPVVsLayer2nd->SetBinContent(bin2D, parameters[1]);
1539       hspectraLGLMPVVsLayer2nd->SetBinError(bin2D, parErrAndRes[1]);
1540       hspectraLGLSigmaVsLayer2nd->SetBinContent(bin2D, parameters[0]);
1541       hspectraLGLSigmaVsLayer2nd->SetBinError(bin2D, parErrAndRes[0]);
1542       hspectraLGGSigmaVsLayer2nd->SetBinContent(bin2D, parameters[3]);
1543       hspectraLGGSigmaVsLayer2nd->SetBinError(bin2D, parErrAndRes[3]);
1544       hMaxLG2nd->Fill(parameters[4]);
1545     }
1546   }
1547   if ( debug > 0)
1548     std::cout << "============================== done fitting 2nd iteration" << std::endl;
1549   int actCh2nd = 0;
1550   double averageScaleUpdated    = calib.GetAverageScaleHigh(actCh2nd);
1551   double meanFWHMHGUpdated      = calib.GetAverageScaleWidthHigh();
1552   double averageScaleLowUpdated = calib.GetAverageScaleLow();
1553   double meanFWHMLGUpdated      = calib.GetAverageScaleWidthLow();
1554   std::cout << "average 1st iteration HG mip: " << averageScale << "\t act. channels: " <<   actCh1st << std::endl;
1555   std::cout << "average 2nd iteration  HG mip: " << averageScaleUpdated << "\t LG mip: " << averageScaleLowUpdated << "\t act. channels: " <<   actCh2nd << std::endl;
1556
1557   RootOutput->cd();
1558
1559   if (IsCalibSaveToFile()){
1560     TString fileCalibPrint = RootOutputName;
1561     fileCalibPrint         = fileCalibPrint.ReplaceAll(".root","_calib.txt");
1562     calib.PrintCalibToFile(fileCalibPrint);
1563   }
1564
1565   TcalibOut->Fill();
1566   TcalibOut->Write();
1567   RootOutput->Close();
1568
1569
1570   RootOutputHist->cd("IndividualCells");
1571     for(ithSpectra=hSpectra.begin(); ithSpectra!=hSpectra.end(); ++ithSpectra){
1572       ithSpectra->second.Write(true);
1573     }
1574   RootOutputHist->cd("IndividualCellsTrigg");
1575     for(ithSpectra=hSpectraTrigg.begin(); ithSpectra!=hSpectraTrigg.end(); ++ithSpectra){
1576       ithSpectra->second.Write(true);
1577     }
1578   RootOutputHist->cd();
1579     hMeanPedHGvsCellID->Write();
1580     hMeanPedHG->Write();
1581     hMeanPedLGvsCellID->Write();
1582     hMeanPedLG->Write();
1583
1584     hspectraHGMeanVsLayer->Write();
1585     hspectraLGMeanVsLayer->Write();
1586     hspectraHGSigmaVsLayer->Write();
1587     hspectraLGSigmaVsLayer->Write();
1588     hspectraHGMaxVsLayer1st->Write();
1589     hspectraHGFWHMVsLayer1st->Write();
1590     hspectraHGLMPVVsLayer1st->Write();
1591     hspectraHGLSigmaVsLayer1st->Write();
1592     hspectraHGGSigmaVsLayer1st->Write();
1593     hMaxHG1st->Write();
1594
1595     hLGHGCorrVsLayer->Write();
1596     hHGLGCorrVsLayer->Write();
1597     hLGHGCorr->Write();
1598     hHGLGCorr->Write();
1599
1600     hspectraHGMaxVsLayer2nd->Write();
1601     hspectraHGFWHMVsLayer2nd->Write();
1602     hspectraHGLMPVVsLayer2nd->Write();
1603     hspectraHGLSigmaVsLayer2nd->Write();
1604     hspectraHGGSigmaVsLayer2nd->Write();
1605     hMaxHG2nd->Write();
1606
1607     hspectraLGMaxVsLayer2nd->Write();
1608     hspectraLGFWHMVsLayer2nd->Write();
1609     hspectraLGLMPVVsLayer2nd->Write();
1610     hspectraLGLSigmaVsLayer2nd->Write();
1611     hspectraLGGSigmaVsLayer2nd->Write();
1612     hMaxLG2nd->Write();
1613
1614     for (int c = 0; c< maxChannelPerLayer; c++)
1615       hHGTileSum[c]->Write();
1616     hmipTriggers->Write();
1617     hSuppresionNoise->Write();
1618     hSuppresionSignal->Write();
1619   // fill calib tree & write it
1620   // close open root files
1621   RootOutputHist->Write();
1622   RootOutputHist->Close();
1623
1624   RootInput->Close();
1625
1626   // Get run info object
1627   std::map<int,RunInfo>::iterator it=ri.find(runNr);
1628
1629   // create directory for plot output
1630   TString outputDirPlots = GetPlotOutputDir();
1631   gSystem->Exec("mkdir -p "+outputDirPlots);
1632
1633   //**********************************************************************
1634   // Create canvases for channel overview plotting
1635   //**********************************************************************
1636   Double_t textSizeRel = 0.035;
1637   StyleSettingsBasics("pdf");
1638   SetPlotStyle();
1639
1640   TCanvas* canvas2DCorr = new TCanvas("canvasCorrPlots","",0,0,1450,1300);  // gives the page size
1641   DefaultCancasSettings( canvas2DCorr, 0.08, 0.13, 0.045, 0.07);
1642
1643   canvas2DCorr->SetLogz(0);
1644   PlotSimple2D( canvas2DCorr, hspectraHGMeanVsLayer, -10000, -10000, textSizeRel, Form("%s/HG_NoiseMean.pdf", outputDirPlots.Data()), it->second, 1, kFALSE, "colz", true);
1645   PlotSimple2D( canvas2DCorr, hspectraHGSigmaVsLayer,-10000, -10000, textSizeRel, Form("%s/HG_NoiseSigma.pdf", outputDirPlots.Data()), it->second, 1,  kFALSE, "colz", true);
1646   PlotSimple2D( canvas2DCorr, hspectraHGMaxVsLayer1st, -10000, -10000, textSizeRel, Form("%s/HG_MaxMip_1st.pdf", outputDirPlots.Data()), it->second, 1, kFALSE, "colz", true, Form( "#LT Max_{HG} #GT = %.1f", averageScale));
1647   PlotSimple2D( canvas2DCorr, hspectraHGFWHMVsLayer1st, -10000, -10000, textSizeRel, Form("%s/HG_FWHMMip_1st.pdf", outputDirPlots.Data()), it->second, 1, kFALSE, "colz", true, Form( "#LT FWHM_{HG} #GT = %.1f", meanFWHMHG));
1648   PlotSimple2D( canvas2DCorr, hspectraHGLMPVVsLayer1st, -10000, -10000, textSizeRel, Form("%s/HG_LandMPVMip_1st.pdf", outputDirPlots.Data()), it->second, 1, kFALSE, "colz", true);
1649   PlotSimple2D( canvas2DCorr, hspectraHGLSigmaVsLayer1st, -10000, -10000, textSizeRel, Form("%s/HG_LandSigMip_1st.pdf", outputDirPlots.Data()), it->second, 1, kFALSE, "colz", true);
1650   PlotSimple2D( canvas2DCorr, hspectraHGGSigmaVsLayer1st, -10000, -10000, textSizeRel, Form("%s/HG_GaussSigMip_1st.pdf", outputDirPlots.Data()), it->second, 1, kFALSE, "colz", true);
1651   PlotSimple2D( canvas2DCorr, hspectraHGMaxVsLayer2nd, -10000, -10000, textSizeRel, Form("%s/HG_MaxMip_2nd.pdf", outputDirPlots.Data()), it->second, 1, kFALSE, "colz", true, Form( "#LT Max_{HG} #GT = %.1f", averageScaleUpdated));
1652   PlotSimple2D( canvas2DCorr, hspectraHGFWHMVsLayer2nd, -10000, -10000, textSizeRel, Form("%s/HG_FWHMMip_2nd.pdf", outputDirPlots.Data()), it->second, 1, kFALSE, "colz", true, Form( "#LT FWHM_{HG} #GT = %.1f", meanFWHMHGUpdated));
1653   PlotSimple2D( canvas2DCorr, hspectraHGLMPVVsLayer2nd, -10000, -10000, textSizeRel, Form("%s/HG_LandMPVMip_2nd.pdf", outputDirPlots.Data()), it->second, 1, kFALSE, "colz", true);
1654   PlotSimple2D( canvas2DCorr, hspectraHGLSigmaVsLayer2nd, -10000, -10000, textSizeRel, Form("%s/HG_LandSigMip_2nd.pdf", outputDirPlots.Data()), it->second, 1, kFALSE, "colz", true);
1655   PlotSimple2D( canvas2DCorr, hspectraHGGSigmaVsLayer2nd, -10000, -10000, textSizeRel, Form("%s/HG_GaussSigMip_2nd.pdf", outputDirPlots.Data()), it->second, 1, kFALSE, "colz", true);
1656   canvas2DCorr->SetLogz(1);
1657   TString drawOpt = "colz";
1658   if (yearData == 2023)
1659     drawOpt = "colz,text";
1660   PlotSimple2D( canvas2DCorr, hmipTriggers, -10000, -10000, textSizeRel, Form("%s/MuonTriggers.pdf", outputDirPlots.Data()), it->second, 1, kFALSE, drawOpt, true, Form( "evt. collected = %d", evts));
1661   PlotSimple2D( canvas2DCorr, hSuppresionNoise, -10000, -10000, textSizeRel, Form("%s/SuppressionNoise.pdf", outputDirPlots.Data()), it->second, 1, kFALSE, drawOpt, true);
1662   PlotSimple2D( canvas2DCorr, hSuppresionSignal, -10000, -10000, textSizeRel, Form("%s/SuppressionSignal.pdf", outputDirPlots.Data()), it->second, 1, kFALSE, drawOpt, true);
1663
1664   canvas2DCorr->SetLogz(0);
1665   PlotSimple2D( canvas2DCorr, hspectraLGMeanVsLayer, -10000, -10000, textSizeRel, Form("%s/LG_NoiseMean.pdf", outputDirPlots.Data()), it->second, 1, kFALSE, "colz", true);
1666   PlotSimple2D( canvas2DCorr, hspectraLGSigmaVsLayer, -10000, -10000, textSizeRel, Form("%s/LG_NoiseSigma.pdf", outputDirPlots.Data()), it->second, 1, kFALSE, "colz", true);
1667   PlotSimple2D( canvas2DCorr, hspectraLGMaxVsLayer2nd, -10000, -10000, textSizeRel, Form("%s/LG_MaxMip_2nd.pdf", outputDirPlots.Data()), it->second, 1, kFALSE, "colz", true, Form( "#LT Max_{LG} #GT = %.1f", averageScaleLowUpdated));
1668   PlotSimple2D( canvas2DCorr, hspectraLGFWHMVsLayer2nd, -10000, -10000, textSizeRel, Form("%s/LG_FWHMMip_2nd.pdf", outputDirPlots.Data()), it->second, 1, kFALSE, "colz", true, Form( "#LT FWHM_{LG} #GT = %.1f", meanFWHMLGUpdated));
1669   PlotSimple2D( canvas2DCorr, hspectraLGLMPVVsLayer2nd, -10000, -10000, textSizeRel, Form("%s/LG_LandMPVMip_2nd.pdf", outputDirPlots.Data()), it->second, 1, kFALSE, "colz", true);
1670   PlotSimple2D( canvas2DCorr, hspectraLGLSigmaVsLayer2nd, -10000, -10000, textSizeRel, Form("%s/LG_LandSigMip_2nd.pdf", outputDirPlots.Data()), it->second, 1, kFALSE, "colz", true);
1671   PlotSimple2D( canvas2DCorr, hspectraLGGSigmaVsLayer2nd, -10000, -10000, textSizeRel, Form("%s/LG_GaussSigMip_2nd.pdf", outputDirPlots.Data()), it->second, 1, kFALSE, "colz", true);
1672
1673   PlotSimple2D( canvas2DCorr, hLGHGCorrVsLayer, -10000, -10000, textSizeRel, Form("%s/LG_HG_Corr.pdf", outputDirPlots.Data()), it->second, 1, kTRUE, "colz", true, Form( "#LT a_{LGHG} #GT = %.1f", avLGHGCorr));
1674   PlotSimple2D( canvas2DCorr, hHGLGCorrVsLayer, -10000, -10000, textSizeRel, Form("%s/HG_LG_Corr.pdf", outputDirPlots.Data()), it->second, 1, kTRUE, "colz", true, Form( "#LT a_{HGLG} #GT = %.1f", avHGLGCorr));
1675
1676   if (ExtPlot > 0){
1677     //***********************************************************************************************************
1678     //********************************* 8 Panel overview plot  **************************************************
1679     //***********************************************************************************************************
1680     //*****************************************************************
1681       // Test beam geometry (beam coming from viewer)
1682       //===========================================================
1683       //||    8 (4)    ||    7 (5)   ||    6 (6)   ||    5 (7)   ||  row 0
1684       //===========================================================
1685       //||    1 (0)    ||    2 (1)   ||    3 (2)   ||    4 (3)   ||  row 1
1686       //===========================================================
1687       //    col 0     col 1       col 2     col  3
1688       // rebuild pad geom in similar way (numbering -1)
1689     //*****************************************************************
1690     TCanvas* canvas8Panel;
1691     TPad* pad8Panel[8];
1692     Double_t topRCornerX[8];
1693     Double_t topRCornerY[8];
1694     Int_t textSizePixel = 30;
1695     Double_t relSize8P[8];
1696     CreateCanvasAndPadsFor8PannelTBPlot(canvas8Panel, pad8Panel,  topRCornerX, topRCornerY, relSize8P, textSizePixel);
1697
1698     TCanvas* canvas8PanelProf;
1699     TPad* pad8PanelProf[8];
1700     Double_t topRCornerXProf[8];
1701     Double_t topRCornerYProf[8];
1702     Double_t relSize8PProf[8];
1703     CreateCanvasAndPadsFor8PannelTBPlot(canvas8PanelProf, pad8PanelProf,  topRCornerXProf, topRCornerYProf, relSize8PProf, textSizePixel, 0.045, "Prof", false);
1704
1705     calib.PrintGlobalInfo();
1706     Double_t maxHG = ReturnMipPlotRangeDepVov(calib.GetVov(),true);
1707     Double_t maxLG = ReturnMipPlotRangeDepVov(calib.GetVov(),false);
1708     std::cout << "plotting single layers" << std::endl;
1709
1710     for (Int_t l = 0; l < setup->GetNMaxLayer()+1; l++){
1711       if (l%10 == 0 && l > 0 && debug > 0)
1712         std::cout << "============================== layer " <<  l << " / " << setup->GetNMaxLayer() << " layers" << std::endl;
1713       if (ExtPlot > 0){
1714         PlotMipWithFitsFullLayer (canvas8Panel,pad8Panel, topRCornerX, topRCornerY, relSize8P, textSizePixel,
1715                                   hSpectra, hSpectraTrigg, setup, true, -100, maxHG, 1.2, l, 0,
1716                                   Form("%s/MIP_HG_Layer%02d.pdf" ,outputDirPlots.Data(), l), it->second);
1717         PlotTriggerPrimWithFitsFullLayer (canvas8Panel,pad8Panel, topRCornerX, topRCornerY, relSize8P, textSizePixel,
1718                                           hSpectraTrigg, setup, averageScale, factorMinTrigg, factorMaxTrigg,
1719                                           0, maxHG*2, 1.2, l, 0, Form("%s/TriggPrimitive_Layer%02d.pdf" ,outputDirPlots.Data(), l), it->second);
1720         PlotCorrWithFitsFullLayer(canvas8PanelProf,pad8PanelProf, topRCornerXProf, topRCornerYProf, relSize8PProf, textSizePixel,
1721                                   hSpectra, setup, false, -20, 800, 1.2, l, 0,
1722                                   Form("%s/LGHG_Corr_Layer%02d.pdf" ,outputDirPlots.Data(), l), it->second);
1723       }
1724       if (ExtPlot > 1){
1725         PlotMipWithFitsFullLayer (canvas8Panel,pad8Panel, topRCornerX, topRCornerY, relSize8P, textSizePixel,
1726                                   hSpectra, hSpectraTrigg, setup, false, -30, maxLG, 1.2, l, 0,
1727                                   Form("%s/MIP_LG_Layer%02d.pdf" ,outputDirPlots.Data(), l), it->second);
1728         PlotCorrWithFitsFullLayer(canvas8PanelProf,pad8PanelProf, topRCornerXProf, topRCornerYProf, relSize8PProf, textSizePixel,
1729                                   hSpectra, setup, true, -100, 4000, 1.2, l, 0,
1730                                   Form("%s/HGLG_Corr_Layer%02d.pdf" ,outputDirPlots.Data(), l), it->second);
1731       }
1732     }
1733     std::cout << "done plotting single layers" << std::endl;
1734   }
1735   return true;
1736 }
1737
1738 //***********************************************************************************************
1739 //*********************** improved scaling calculation function *********************************
1740 //***********************************************************************************************
1741 bool Analyses::GetImprovedScaling(void){
1742   std::cout<<"GetImprovedScaling"<<std::endl;
1743
1744   std::map<int,RunInfo> ri=readRunInfosFromFile(RunListInputName.Data(),debug,0);
1745
1746   std::map<int,TileSpectra> hSpectra;
1747   std::map<int,TileSpectra> hSpectraTrigg;
1748   std::map<int, TileSpectra>::iterator ithSpectra;
1749   std::map<int, TileSpectra>::iterator ithSpectraTrigg;
1750
1751   std::cout << "Additional Output with histos being created: " << RootOutputNameHist.Data() << std::endl;
1752   if(Overwrite){
1753     std::cout << "recreating file with hists" << std::endl;
1754     RootOutputHist = new TFile(RootOutputNameHist.Data(),"RECREATE");
1755   } else{
1756     std::cout << "newly creating file with hists" << std::endl;
1757     RootOutputHist = new TFile(RootOutputNameHist.Data(),"CREATE");
1758   }
1759
1760   ROOT::Math::MinimizerOptions::SetDefaultMinimizer("Minuit2", "Migrad");
1761   // setup trigger sel
1762   double factorMinTrigg   = 0.8;
1763   double factorMaxTrigg   = 2.5;
1764   if (yearData == 2023){
1765     factorMinTrigg    = 0.9;
1766     factorMaxTrigg    = 2.;
1767   }
1768
1769   RootOutputHist->mkdir("IndividualCells");
1770   RootOutputHist->mkdir("IndividualCellsTrigg");
1771   RootOutputHist->cd("IndividualCellsTrigg");
1772   //***********************************************************************************************
1773   //************************* first pass over tree to extract spectra *****************************
1774   //***********************************************************************************************
1775   TcalibIn->GetEntry(0);
1776   int evts=TdataIn->GetEntries();
1777   int runNr = -1;
1778   int actChI  = 0;
1779   double averageScale     = calib.GetAverageScaleHigh(actChI);
1780   double averageScaleLow  = calib.GetAverageScaleLow();
1781   std::cout << "average HG mip: " << averageScale << " LG mip: "<< averageScaleLow << "\t act. ch: "<< actChI << std::endl;
1782   for(int i=0; i<evts; i++){
1783     TdataIn->GetEntry(i);
1784     if (i == 0)runNr = event.GetRunNumber();
1785     TdataIn->GetEntry(i);
1786     if (i%5000 == 0 && i > 0 && debug > 0) std::cout << "Reading " <<  i << " / " << evts << " events" << std::endl;
1787     for(int j=0; j<event.GetNTiles(); j++){
1788       Caen* aTile=(Caen*)event.GetTile(j);
1789       if (i == 0 && debug > 2) std::cout << ((TString)setup->DecodeCellID(aTile->GetCellID())).Data() << std::endl;
1790       long currCellID = aTile->GetCellID();
1791
1792       // read current tile
1793       ithSpectraTrigg=hSpectraTrigg.find(aTile->GetCellID());
1794       double hgCorr = aTile->GetADCHigh()-calib.GetPedestalMeanH(aTile->GetCellID());
1795       double lgCorr = aTile->GetADCLow()-calib.GetPedestalMeanL(aTile->GetCellID());
1796
1797       // estimate local muon trigger
1798       bool localMuonTrigg = event.InspectIfLocalMuonTrigg(currCellID, averageScale, factorMinTrigg, factorMaxTrigg);
1799
1800       if(ithSpectraTrigg!=hSpectraTrigg.end()){
1801         ithSpectraTrigg->second.FillTrigger(aTile->GetLocalTriggerPrimitive());
1802       } else {
1803         RootOutputHist->cd("IndividualCellsTrigg");
1804         hSpectraTrigg[currCellID]=TileSpectra("mipTrigg",currCellID,calib.GetTileCalib(currCellID),debug);
1805         hSpectraTrigg[currCellID].FillTrigger(aTile->GetLocalTriggerPrimitive());;
1806         RootOutput->cd();
1807       }
1808
1809       ithSpectra=hSpectra.find(aTile->GetCellID());
1810       if (ithSpectra!=hSpectra.end()){
1811         ithSpectra->second.FillSpectra(lgCorr,hgCorr);
1812         if (hgCorr > 3*calib.GetPedestalSigH(currCellID) && lgCorr > 3*calib.GetPedestalSigL(currCellID) && hgCorr < 3900 )
1813           ithSpectra->second.FillCorr(lgCorr,hgCorr);
1814       } else {
1815         RootOutputHist->cd("IndividualCells");
1816         hSpectra[currCellID]=TileSpectra("mip1st",currCellID,calib.GetTileCalib(currCellID),debug);
1817         hSpectra[aTile->GetCellID()].FillSpectra(lgCorr,hgCorr);;
1818         if (hgCorr > 3*calib.GetPedestalSigH(aTile->GetCellID()) && lgCorr > 3*calib.GetPedestalSigL(aTile->GetCellID() && hgCorr < 3900) )
1819           hSpectra[aTile->GetCellID()].FillCorr(lgCorr,hgCorr);;
1820
1821         RootOutput->cd();
1822       }
1823
1824
1825       // only fill tile spectra if 4 surrounding tiles on average are compatible with muon response
1826       if (localMuonTrigg){
1827         aTile->SetLocalTriggerBit(1);
1828         ithSpectraTrigg=hSpectraTrigg.find(aTile->GetCellID());
1829         ithSpectraTrigg->second.FillSpectra(lgCorr,hgCorr);
1830         if (hgCorr > 3*calib.GetPedestalSigH(currCellID) && lgCorr > 3*calib.GetPedestalSigL(currCellID) && hgCorr < 3900 )
1831           ithSpectraTrigg->second.FillCorr(lgCorr,hgCorr);
1832       }
1833     }
1834     RootOutput->cd();
1835     TdataOut->Fill();
1836   }
1837   TdataOut->Write();
1838   TsetupIn->CloneTree()->Write();
1839
1840   //***********************************************************************************************
1841   //***** Monitoring histos for fits results of 2nd iteration ******************
1842   //***********************************************************************************************
1843   RootOutputHist->cd();
1844   int maxChannelPerLayer        = (setup->GetNMaxColumn()+1)*(setup->GetNMaxRow()+1);
1845   // monitoring trigger
1846   TH2D* hmipTriggers              = new TH2D( "hmipTriggers","muon triggers; layer; brd channel; counts ",
1847                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
1848   hmipTriggers->SetDirectory(0);
1849   TH2D* hSuppresionNoise          = new TH2D( "hSuppresionNoise","S/B noise region; layer; brd channel; S/B noise region",
1850                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
1851   hSuppresionNoise->SetDirectory(0);
1852   TH2D* hSuppresionSignal         = new TH2D( "hSuppresionSignal","S/B signal region; layer; brd channel; S/B signal region",
1853                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
1854   hSuppresionSignal->SetDirectory(0);
1855
1856   // monitoring 2nd iteration mip fits
1857   TH2D* hspectraHGMaxVsLayer   = new TH2D( "hspectraHGMaxVsLayer","Max High Gain; layer; brd channel; Max_{HG, 2^{nd}} (arb. units) ",
1858                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
1859   hspectraHGMaxVsLayer->SetDirectory(0);
1860   TH2D* hspectraHGFWHMVsLayer   = new TH2D( "hspectraHGFWHMVsLayer","FWHM High Gain; layer; brd channel; FWHM_{HG, 2^{nd}} (arb. units) ",
1861                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
1862   hspectraHGFWHMVsLayer->SetDirectory(0);
1863   TH2D* hspectraHGLMPVVsLayer   = new TH2D( "hspectraHGLMPVVsLayer","MPV High Gain; layer; brd channel; MPV_{HG, 2^{nd}} (arb. units) ",
1864                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
1865   hspectraHGLMPVVsLayer->SetDirectory(0);
1866   TH2D* hspectraHGLSigmaVsLayer = new TH2D( "hspectraHGLSigmaVsLayer","Sigma Landau High Gain; layer; brd channel; #sigma_{L,HG, 2^{nd}} (arb. units) ",
1867                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
1868   hspectraHGLSigmaVsLayer->SetDirectory(0);
1869   TH2D* hspectraHGGSigmaVsLayer = new TH2D( "hspectraHGGSigmaVsLayer","Sigma Gauss High Gain; layer; brd channel; #sigma_{G,HG, 2^{nd}} (arb. units) ",
1870                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
1871   hspectraHGGSigmaVsLayer->SetDirectory(0);
1872   TH2D* hspectraLGMaxVsLayer   = new TH2D( "hspectraLGMaxVsLayer","Max High Gain; layer; brd channel; Max_{LG, 2^{nd}} (arb. units) ",
1873                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
1874   hspectraLGMaxVsLayer->SetDirectory(0);
1875   TH2D* hspectraLGFWHMVsLayer   = new TH2D( "hspectraLGFWHMVsLayer","FWHM High Gain; layer; brd channel; FWHM_{LG, 2^{nd}} (arb. units) ",
1876                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
1877   hspectraLGFWHMVsLayer->SetDirectory(0);
1878   TH2D* hspectraLGLMPVVsLayer   = new TH2D( "hspectraLGLMPVVsLayer","MPV High Gain; layer; brd channel; MPV_{LG, 2^{nd}} (arb. units) ",
1879                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
1880   hspectraLGLMPVVsLayer->SetDirectory(0);
1881   TH2D* hspectraLGLSigmaVsLayer = new TH2D( "hspectraLGLSigmaVsLayer","Sigma Landau High Gain; layer; brd channel; #sigma_{L,LG, 2^{nd}} (arb. units) ",
1882                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
1883   hspectraLGLSigmaVsLayer->SetDirectory(0);
1884   TH2D* hspectraLGGSigmaVsLayer = new TH2D( "hspectraLGGSigmaVsLayer","Sigma Gauss High Gain; layer; brd channel; #sigma_{G,LG, 2^{nd}} (arb. units) ",
1885                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
1886   hspectraLGGSigmaVsLayer->SetDirectory(0);
1887
1888   TH1D* hMaxHG             = new TH1D( "hMaxHG","Max High Gain ;Max_{HG} (arb. units) ; counts ",
1889                                             2000, -0.5, 2000-0.5);
1890   hMaxHG->SetDirectory(0);
1891   TH1D* hMaxLG             = new TH1D( "hMaxLG","Max Low Gain ;Max_{LG} (arb. units) ; counts ",
1892                                             400, -0.5, 400-0.5);
1893   hMaxLG->SetDirectory(0);
1894
1895
1896   int currCells = 0;
1897   double* parameters    = new double[6];
1898   double* parErrAndRes  = new double[6];
1899   bool isGood;
1900   double meanSB_NoiseR  = 0;
1901   double meanSB_SigR    = 0;
1902   if ( debug > 0)
1903     std::cout << "============================== start fitting improved iteration" << std::endl;
1904
1905   for(ithSpectraTrigg=hSpectraTrigg.begin(); ithSpectraTrigg!=hSpectraTrigg.end(); ++ithSpectraTrigg){
1906     if (currCells%20 == 0 && currCells > 0 && debug > 0)
1907       std::cout << "============================== cell " <<  currCells << " / " << hSpectraTrigg.size() << " cells" << std::endl;
1908     currCells++;
1909     for (int p = 0; p < 6; p++){
1910       parameters[p]   = 0;
1911       parErrAndRes[p] = 0;
1912     }
1913     isGood=ithSpectraTrigg->second.FitMipHG(parameters, parErrAndRes, debug, yearData, true, calib.GetVov(), averageScale);
1914
1915     long cellID     = ithSpectraTrigg->second.GetCellID();
1916     int layer       = setup->GetLayer(cellID);
1917     int chInLayer   = setup->GetChannelInLayer(cellID);
1918     int bin2D       = hspectraHGMaxVsLayer->FindBin(layer,chInLayer);
1919
1920     Int_t binNLow   = ithSpectraTrigg->second.GetHG()->FindBin(-1*calib.GetPedestalSigH(cellID));
1921     Int_t binNHigh  = ithSpectraTrigg->second.GetHG()->FindBin(3*calib.GetPedestalSigH(cellID));
1922     Int_t binSHigh  = ithSpectraTrigg->second.GetHG()->FindBin(3800);
1923
1924     double S_NoiseR = ithSpectraTrigg->second.GetHG()->Integral(binNLow, binNHigh);
1925     double S_SigR   = ithSpectraTrigg->second.GetHG()->Integral(binNHigh, binSHigh);
1926
1927     ithSpectra      = hSpectra.find(cellID);
1928     double B_NoiseR = ithSpectra->second.GetHG()->Integral(binNLow , binNHigh);
1929     double B_SigR   = ithSpectra->second.GetHG()->Integral(binNHigh, binSHigh);
1930
1931     double SB_NoiseR  = (B_NoiseR != 0.) ? S_NoiseR/B_NoiseR : 0;
1932     double SB_SigR    = (B_SigR != 0.) ? S_SigR/B_SigR : 0;
1933
1934     meanSB_NoiseR += SB_NoiseR;
1935     meanSB_SigR += SB_SigR;
1936
1937     hmipTriggers->SetBinContent(bin2D, ithSpectraTrigg->second.GetHG()->GetEntries());
1938     hSuppresionNoise->SetBinContent(bin2D, SB_NoiseR);
1939     hSuppresionSignal->SetBinContent(bin2D, SB_SigR);
1940     if (isGood){
1941       hspectraHGMaxVsLayer->SetBinContent(bin2D, parameters[4]);
1942       hspectraHGFWHMVsLayer->SetBinContent(bin2D, parameters[5]);
1943       hspectraHGLMPVVsLayer->SetBinContent(bin2D, parameters[1]);
1944       hspectraHGLMPVVsLayer->SetBinError(bin2D, parErrAndRes[1]);
1945       hspectraHGLSigmaVsLayer->SetBinContent(bin2D, parameters[0]);
1946       hspectraHGLSigmaVsLayer->SetBinError(bin2D, parErrAndRes[0]);
1947       hspectraHGGSigmaVsLayer->SetBinContent(bin2D, parameters[3]);
1948       hspectraHGGSigmaVsLayer->SetBinError(bin2D, parErrAndRes[3]);
1949       hMaxHG->Fill(parameters[4]);
1950     }
1951     for (int p = 0; p < 6; p++){
1952       parameters[p]   = 0;
1953       parErrAndRes[p] = 0;
1954     }
1955     isGood=ithSpectraTrigg->second.FitMipLG(parameters, parErrAndRes, debug, yearData, true, averageScaleLow);
1956     if (isGood){
1957       hspectraLGMaxVsLayer->SetBinContent(bin2D, parameters[4]);
1958       hspectraLGFWHMVsLayer->SetBinContent(bin2D, parameters[5]);
1959       hspectraLGLMPVVsLayer->SetBinContent(bin2D, parameters[1]);
1960       hspectraLGLMPVVsLayer->SetBinError(bin2D, parErrAndRes[1]);
1961       hspectraLGLSigmaVsLayer->SetBinContent(bin2D, parameters[0]);
1962       hspectraLGLSigmaVsLayer->SetBinError(bin2D, parErrAndRes[0]);
1963       hspectraLGGSigmaVsLayer->SetBinContent(bin2D, parameters[3]);
1964       hspectraLGGSigmaVsLayer->SetBinError(bin2D, parErrAndRes[3]);
1965       hMaxLG->Fill(parameters[4]);
1966     }
1967   }
1968   if ( debug > 0)
1969     std::cout << "============================== done fitting improved iteration" << std::endl;
1970
1971
1972   meanSB_NoiseR = meanSB_NoiseR/hSpectraTrigg.size();
1973   meanSB_SigR   = meanSB_SigR/hSpectraTrigg.size();
1974
1975   RootOutput->cd();
1976
1977   if (IsCalibSaveToFile()){
1978     TString fileCalibPrint = RootOutputName;
1979     fileCalibPrint         = fileCalibPrint.ReplaceAll(".root","_calib.txt");
1980     calib.PrintCalibToFile(fileCalibPrint);
1981   }
1982
1983   TcalibOut->Fill();
1984   TcalibOut->Write();
1985   int actChA                     = 0;
1986   double averageScaleUpdated     = calib.GetAverageScaleHigh(actChA);
1987   double averageScaleUpdatedLow  = calib.GetAverageScaleLow();
1988   double meanFWHMHG              = calib.GetAverageScaleWidthHigh();
1989   double meanFWHMLG              = calib.GetAverageScaleWidthLow();
1990
1991   std::cout << "average input HG mip: " << averageScale << " LG mip: "<< averageScaleLow << "\t act. ch: "<< actChI<< std::endl;
1992   std::cout << "average updated HG mip: " << averageScaleUpdated << " LG mip: "<< averageScaleUpdatedLow << "\t act. ch: "<< actChA<< std::endl;
1993
1994   RootOutput->Close();
1995
1996
1997   RootOutputHist->cd("IndividualCellsTrigg");
1998     for(ithSpectra=hSpectraTrigg.begin(); ithSpectra!=hSpectraTrigg.end(); ++ithSpectra){
1999       ithSpectra->second.Write(true);
2000     }
2001   RootOutputHist->cd();
2002
2003     hspectraHGMaxVsLayer->Write();
2004     hspectraHGFWHMVsLayer->Write();
2005     hspectraHGLMPVVsLayer->Write();
2006     hspectraHGLSigmaVsLayer->Write();
2007     hspectraHGGSigmaVsLayer->Write();
2008     hMaxHG->Write();
2009
2010     hspectraLGMaxVsLayer->Write();
2011     hspectraLGFWHMVsLayer->Write();
2012     hspectraLGLMPVVsLayer->Write();
2013     hspectraLGLSigmaVsLayer->Write();
2014     hspectraLGGSigmaVsLayer->Write();
2015     hMaxLG->Write();
2016
2017     hmipTriggers->Write();
2018     hSuppresionNoise->Write();
2019     hSuppresionSignal->Write();
2020   // fill calib tree & write it
2021   // close open root files
2022   RootOutputHist->Write();
2023   RootOutputHist->Close();
2024
2025   RootInput->Close();
2026
2027   // Get run info object
2028   std::map<int,RunInfo>::iterator it=ri.find(runNr);
2029
2030   // create directory for plot output
2031   TString outputDirPlots = GetPlotOutputDir();
2032   gSystem->Exec("mkdir -p "+outputDirPlots);
2033
2034   //**********************************************************************
2035   // Create canvases for channel overview plotting
2036   //**********************************************************************
2037   Double_t textSizeRel = 0.035;
2038   StyleSettingsBasics("pdf");
2039   SetPlotStyle();
2040
2041   TCanvas* canvas2DCorr = new TCanvas("canvasCorrPlots","",0,0,1450,1300);  // gives the page size
2042   DefaultCancasSettings( canvas2DCorr, 0.08, 0.13, 0.045, 0.07);
2043
2044   canvas2DCorr->SetLogz(0);
2045   PlotSimple2D( canvas2DCorr, hspectraHGMaxVsLayer, -10000, -10000, textSizeRel, Form("%s/HG_MaxMip.pdf", outputDirPlots.Data()), it->second, 1, kFALSE, "colz", true, Form( "#LT Max_{HG} #GT = %.1f", averageScaleUpdated) );
2046   PlotSimple2D( canvas2DCorr, hspectraHGFWHMVsLayer, -10000, -10000, textSizeRel, Form("%s/HG_FWHMMip.pdf", outputDirPlots.Data()), it->second, 1, kFALSE, "colz", true, Form( "#LT FWHM_{HG} #GT = %.1f", meanFWHMHG));
2047   PlotSimple2D( canvas2DCorr, hspectraHGLMPVVsLayer, -10000, -10000, textSizeRel, Form("%s/HG_LandMPVMip.pdf", outputDirPlots.Data()), it->second, 1, kFALSE, "colz", true);
2048   PlotSimple2D( canvas2DCorr, hspectraHGLSigmaVsLayer, -10000, -10000, textSizeRel, Form("%s/HG_LandSigMip.pdf", outputDirPlots.Data()), it->second, 1, kFALSE, "colz", true);
2049   PlotSimple2D( canvas2DCorr, hspectraHGGSigmaVsLayer, -10000, -10000, textSizeRel, Form("%s/HG_GaussSigMip.pdf", outputDirPlots.Data()), it->second, 1, kFALSE, "colz", true);
2050   canvas2DCorr->SetLogz(1);
2051   TString drawOpt = "colz";
2052   if (yearData == 2023)
2053     drawOpt = "colz,text";
2054   PlotSimple2D( canvas2DCorr, hmipTriggers, -10000, -10000, textSizeRel, Form("%s/MuonTriggers.pdf", outputDirPlots.Data()), it->second, 1, kFALSE, drawOpt, true, Form( "evt. collected = %d", evts));
2055   PlotSimple2D( canvas2DCorr, hSuppresionNoise, -10000, -10000, textSizeRel, Form("%s/SuppressionNoise.pdf", outputDirPlots.Data()), it->second, 1, kFALSE, drawOpt, true, Form( "#LT S/B noise #GT = %.3f", meanSB_NoiseR));
2056   PlotSimple2D( canvas2DCorr, hSuppresionSignal, -10000, -10000, textSizeRel, Form("%s/SuppressionSignal.pdf", outputDirPlots.Data()), it->second, 1, kFALSE, drawOpt, true, Form( "#LT S/B signal #GT = %.3f", meanSB_SigR));
2057
2058   canvas2DCorr->SetLogz(0);
2059   PlotSimple2D( canvas2DCorr, hspectraLGMaxVsLayer, -10000, -10000, textSizeRel, Form("%s/LG_MaxMip.pdf", outputDirPlots.Data()), it->second, 1, kFALSE, "colz", true, Form( "#LT Max_{LG} #GT = %.1f", averageScaleUpdatedLow));
2060   PlotSimple2D( canvas2DCorr, hspectraLGFWHMVsLayer, -10000, -10000, textSizeRel, Form("%s/LG_FWHMMip.pdf", outputDirPlots.Data()), it->second, 1, kFALSE, "colz", true, Form( "#LT FWHM_{LG} #GT = %.1f", meanFWHMLG));
2061   PlotSimple2D( canvas2DCorr, hspectraLGLMPVVsLayer, -10000, -10000, textSizeRel, Form("%s/LG_LandMPVMip.pdf", outputDirPlots.Data()), it->second, 1, kFALSE, "colz", true);
2062   PlotSimple2D( canvas2DCorr, hspectraLGLSigmaVsLayer, -10000, -10000, textSizeRel, Form("%s/LG_LandSigMip.pdf", outputDirPlots.Data()), it->second, 1, kFALSE, "colz", true);
2063   PlotSimple2D( canvas2DCorr, hspectraLGGSigmaVsLayer, -10000, -10000, textSizeRel, Form("%s/LG_GaussSigMip.pdf", outputDirPlots.Data()), it->second, 1, kFALSE, "colz", true);
2064
2065   //***********************************************************************************************************
2066   //********************************* 8 Panel overview plot  **************************************************
2067   //***********************************************************************************************************
2068   //*****************************************************************
2069     // Test beam geometry (beam coming from viewer)
2070     //===========================================================
2071     //||    8 (4)    ||    7 (5)   ||    6 (6)   ||    5 (7)   ||  row 0
2072     //===========================================================
2073     //||    1 (0)    ||    2 (1)   ||    3 (2)   ||    4 (3)   ||  row 1
2074     //===========================================================
2075     //    col 0     col 1       col 2     col  3
2076     // rebuild pad geom in similar way (numbering -1)
2077   //*****************************************************************
2078   TCanvas* canvas8Panel;
2079   TPad* pad8Panel[8];
2080   Double_t topRCornerX[8];
2081   Double_t topRCornerY[8];
2082   Int_t textSizePixel = 30;
2083   Double_t relSize8P[8];
2084   CreateCanvasAndPadsFor8PannelTBPlot(canvas8Panel, pad8Panel,  topRCornerX, topRCornerY, relSize8P, textSizePixel);
2085
2086   calib.PrintGlobalInfo();
2087   Double_t maxHG = ReturnMipPlotRangeDepVov(calib.GetVov(),true);
2088   Double_t maxLG = ReturnMipPlotRangeDepVov(calib.GetVov(),false);
2089   std::cout << "plotting single layers" << std::endl;
2090   for (Int_t l = 0; l < setup->GetNMaxLayer()+1; l++){
2091     if (l%10 == 0 && l > 0 && debug > 0)
2092       std::cout << "============================== layer " <<  l << " / " << setup->GetNMaxLayer() << " layers" << std::endl;
2093     PlotMipWithFitsFullLayer (canvas8Panel,pad8Panel, topRCornerX, topRCornerY, relSize8P, textSizePixel,
2094                               hSpectra, hSpectraTrigg, setup, true, -100, maxHG, 1.2, l, 0,
2095                               Form("%s/MIP_HG_Layer%02d.pdf" ,outputDirPlots.Data(), l), it->second);
2096     PlotMipWithFitsFullLayer (canvas8Panel,pad8Panel, topRCornerX, topRCornerY, relSize8P, textSizePixel,
2097                               hSpectra, hSpectraTrigg, setup, false, -20, maxLG, 1.2, l, 0,
2098                               Form("%s/MIP_LG_Layer%02d.pdf" ,outputDirPlots.Data(), l), it->second);
2099     PlotTriggerPrimWithFitsFullLayer (canvas8Panel,pad8Panel, topRCornerX, topRCornerY, relSize8P, textSizePixel,
2100                                       hSpectraTrigg, setup, averageScale, factorMinTrigg, factorMaxTrigg,
2101                                       0, maxHG*2, 1.2, l, 0, Form("%s/TriggPrimitive_Layer%02d.pdf" ,outputDirPlots.Data(), l), it->second);
2102   }
2103   std::cout << "done plotting" << std::endl;
2104
2105
2106   return true;
2107 }
2108
2109
2110
2111 //***********************************************************************************************
2112 //*********************** improved scaling calculation function *********************************
2113 //***********************************************************************************************
2114 bool Analyses::GetNoiseSampleAndRefitPedestal(void){
2115   std::cout<<"GetNoiseSampleAndRefitPedestal"<<std::endl;
2116
2117   std::map<int,RunInfo> ri=readRunInfosFromFile(RunListInputName.Data(),debug,0);
2118
2119   std::map<int,TileSpectra> hSpectra;
2120   std::map<int,TileSpectra> hSpectraTrigg;
2121   std::map<int, TileSpectra>::iterator ithSpectra;
2122   std::map<int, TileSpectra>::iterator ithSpectraTrigg;
2123
2124   std::cout << "Additional Output with histos being created: " << RootOutputNameHist.Data() << std::endl;
2125   if(Overwrite){
2126     std::cout << "recreating file with hists" << std::endl;
2127     RootOutputHist = new TFile(RootOutputNameHist.Data(),"RECREATE");
2128   } else{
2129     std::cout << "newly creating file with hists" << std::endl;
2130     RootOutputHist = new TFile(RootOutputNameHist.Data(),"CREATE");
2131   }
2132
2133   // setup trigger sel
2134   double factorMinTrigg   = 0.5;
2135   if(yearData == 2023)
2136     factorMinTrigg        = 0.1;
2137   // create HG and LG histo's per channel
2138   TH2D* hspectraHGvsCellID      = new TH2D( "hNoiseTriggeredSpectraHG_vsCellID","Noise trigg ADC spectrum High Gain vs CellID; cell ID; ADC_{HG} (arb. units) ",
2139                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5, 4000,0,4000);
2140   hspectraHGvsCellID->SetDirectory(0);
2141   TH2D* hspectraLGvsCellID      = new TH2D( "hNoiseTriggeredSpectraLG_vsCellID","Noise trigg ADC spectrum Low Gain vs CellID; cell ID; ADC_{LG} (arb. units)  ",
2142                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5, 4000,0,4000);
2143   hspectraLGvsCellID->SetDirectory(0);
2144
2145
2146   RootOutputHist->mkdir("IndividualCells");
2147   RootOutputHist->mkdir("IndividualCellsTrigg");
2148   RootOutputHist->cd("IndividualCellsTrigg");
2149
2150   //***********************************************************************************************
2151   //************************* first pass over tree to extract spectra *****************************
2152   //***********************************************************************************************
2153   TcalibIn->GetEntry(0);
2154   int evts=TdataIn->GetEntries();
2155   int runNr = -1;
2156   int actCh = 0;
2157   double averageScale     = calib.GetAverageScaleHigh(actCh);
2158   double averageScaleLow  = calib.GetAverageScaleLow();
2159   std::cout << "average HG mip: " << averageScale << " LG mip: "<< averageScaleLow << std::endl;
2160   for(int i=0; i<evts; i++){
2161     TdataIn->GetEntry(i);
2162     if (i == 0)runNr = event.GetRunNumber();
2163     TdataIn->GetEntry(i);
2164     if (i%5000 == 0 && i > 0 && debug > 0) std::cout << "Reading " <<  i << " / " << evts << " events" << std::endl;
2165     for(int j=0; j<event.GetNTiles(); j++){
2166       Caen* aTile=(Caen*)event.GetTile(j);
2167       if (i == 0 && debug > 2) std::cout << ((TString)setup->DecodeCellID(aTile->GetCellID())).Data() << std::endl;
2168       long currCellID = aTile->GetCellID();
2169
2170       // read current tile
2171       ithSpectraTrigg=hSpectraTrigg.find(aTile->GetCellID());
2172       // estimate local muon trigger
2173       bool localNoiseTrigg = event.InspectIfNoiseTrigg(currCellID, averageScale, factorMinTrigg);
2174
2175       if(ithSpectraTrigg!=hSpectraTrigg.end()){
2176         ithSpectraTrigg->second.FillTrigger(aTile->GetLocalTriggerPrimitive());
2177       } else {
2178         RootOutputHist->cd("IndividualCellsTrigg");
2179         hSpectraTrigg[currCellID]=TileSpectra("mipTrigg",currCellID,calib.GetTileCalib(currCellID),debug);
2180         hSpectraTrigg[currCellID].FillTrigger(aTile->GetLocalTriggerPrimitive());;
2181         RootOutput->cd();
2182       }
2183
2184       ithSpectra=hSpectra.find(aTile->GetCellID());
2185       if (ithSpectra!=hSpectra.end()){
2186         ithSpectra->second.FillSpectra(aTile->GetADCLow(),aTile->GetADCHigh());
2187       } else {
2188         RootOutputHist->cd("IndividualCells");
2189         hSpectra[currCellID]=TileSpectra("mip1st",currCellID,calib.GetTileCalib(currCellID),debug);
2190         hSpectra[aTile->GetCellID()].FillSpectra(aTile->GetADCLow(),aTile->GetADCHigh());;
2191
2192         RootOutput->cd();
2193       }
2194
2195       // only fill tile spectra if X surrounding tiles on average are compatible with pure noise
2196       if (localNoiseTrigg){
2197         aTile->SetLocalTriggerBit(2);
2198         ithSpectraTrigg=hSpectraTrigg.find(aTile->GetCellID());
2199         ithSpectraTrigg->second.FillSpectra(aTile->GetADCLow(),aTile->GetADCHigh());
2200
2201         hspectraHGvsCellID->Fill(aTile->GetCellID(), aTile->GetADCHigh());
2202         hspectraLGvsCellID->Fill(aTile->GetCellID(), aTile->GetADCLow());
2203       }
2204     }
2205     RootOutput->cd();
2206     TdataOut->Fill();
2207   }
2208   TdataOut->Write();
2209   TsetupIn->CloneTree()->Write();
2210
2211   //***********************************************************************************************
2212   //***** Monitoring histos for fits results of 2nd iteration ******************
2213   //***********************************************************************************************
2214   RootOutputHist->cd();
2215   int maxChannelPerLayer        = (setup->GetNMaxColumn()+1)*(setup->GetNMaxRow()+1);
2216   // monitoring trigger
2217   TH2D* hnoiseTriggers              = new TH2D( "hnoiseTriggers","muon triggers; layer; brd channel; counts ",
2218                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
2219   hnoiseTriggers->SetDirectory(0);
2220   TH1D* hMeanPedHGvsCellID      = new TH1D( "hMeanPedHG_vsCellID","mean Ped High Gain vs CellID ; cell ID; #mu_{noise, HG} (arb. units) ",
2221                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5);
2222   hMeanPedHGvsCellID->SetDirectory(0);
2223   TH2D* hspectraHGMeanVsLayer   = new TH2D( "hspectraHGMeanVsLayer","Mean Ped High Gain vs CellID; layer; brd channel; #mu_{Ped HG} (arb. units) ",
2224                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
2225   hspectraHGMeanVsLayer->SetDirectory(0);
2226   TH2D* hspectraHGSigmaVsLayer  = new TH2D( "hspectraHGSigmaVsLayer","Mean Ped High Gain vs CellID; layer; brd channel; #sigma_{Ped HG} (arb. units) ",
2227                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
2228   hspectraHGSigmaVsLayer->SetDirectory(0);
2229   TH1D* hMeanPedLGvsCellID      = new TH1D( "hMeanPedLG_vsCellID","mean Ped Low Gain vs CellID ; cell ID; #mu_{noise, LG} (arb. units) ",
2230                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5);
2231   hMeanPedLGvsCellID->SetDirectory(0);
2232   TH2D* hspectraLGMeanVsLayer   = new TH2D( "hspectraLGMeanVsLayer","Mean Ped Low Gain vs CellID; layer; brd channel; #mu_{PED LG} (arb. units) ",
2233                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
2234   hspectraLGMeanVsLayer->SetDirectory(0);
2235   TH2D* hspectraLGSigmaVsLayer  = new TH2D( "hspectraLGSigmaVsLayer","Mean Ped Low Gain vs CellID; layer; brd channel; #sigma_{Ped LG} (arb. units)",
2236                                             setup->GetNMaxLayer()+1, -0.5, setup->GetNMaxLayer()+1-0.5, maxChannelPerLayer, -0.5, maxChannelPerLayer-0.5);
2237   hspectraLGSigmaVsLayer->SetDirectory(0);
2238
2239   if ( debug > 0)
2240     std::cout << "============================== starting fitting" << std::endl;
2241
2242   int currCells = 0;
2243   double* parameters    = new double[6];
2244   for(ithSpectraTrigg=hSpectraTrigg.begin(); ithSpectraTrigg!=hSpectraTrigg.end(); ++ithSpectraTrigg){
2245     for (int p = 0; p < 6; p++){
2246       parameters[p]   = 0;
2247     }
2248
2249     if (currCells%20 == 0 && currCells > 0 && debug > 0)
2250       std::cout << "============================== cell " <<  currCells << " / " << hSpectraTrigg.size() << " cells" << std::endl;
2251     currCells++;
2252     if ( debug > 2) std::cout << ((TString)setup->DecodeCellID(ithSpectraTrigg->second.GetCellID())).Data() << std::endl;
2253     ithSpectraTrigg->second.FitNoise(parameters, yearData, true);
2254     hMeanPedHGvsCellID->SetBinContent(hMeanPedHGvsCellID->FindBin(ithSpectraTrigg->second.GetCellID()), parameters[4]);
2255     hMeanPedHGvsCellID->SetBinError  (hMeanPedHGvsCellID->FindBin(ithSpectraTrigg->second.GetCellID()), parameters[6]);
2256     hMeanPedLGvsCellID->SetBinContent(hMeanPedLGvsCellID->FindBin(ithSpectraTrigg->second.GetCellID()), parameters[0]);
2257     hMeanPedLGvsCellID->SetBinError  (hMeanPedLGvsCellID->FindBin(ithSpectraTrigg->second.GetCellID()), parameters[2]);
2258
2259     int layer     = setup->GetLayer(ithSpectraTrigg->second.GetCellID());
2260     int chInLayer = setup->GetChannelInLayer(ithSpectraTrigg->second.GetCellID());
2261
2262     hspectraHGMeanVsLayer->SetBinContent(hspectraHGMeanVsLayer->FindBin(layer,chInLayer), parameters[4]);
2263     hspectraHGMeanVsLayer->SetBinError(hspectraHGMeanVsLayer->FindBin(layer,chInLayer), parameters[5]);
2264     hspectraHGSigmaVsLayer->SetBinContent(hspectraHGSigmaVsLayer->FindBin(layer,chInLayer), parameters[6]);
2265     hspectraHGSigmaVsLayer->SetBinError(hspectraHGSigmaVsLayer->FindBin(layer,chInLayer), parameters[7]);
2266     hspectraLGMeanVsLayer->SetBinContent(hspectraLGMeanVsLayer->FindBin(layer,chInLayer), parameters[0]);
2267     hspectraLGMeanVsLayer->SetBinError(hspectraLGMeanVsLayer->FindBin(layer,chInLayer), parameters[1]);
2268     hspectraLGSigmaVsLayer->SetBinContent(hspectraLGSigmaVsLayer->FindBin(layer,chInLayer), parameters[2]);
2269     hspectraLGSigmaVsLayer->SetBinError(hspectraLGSigmaVsLayer->FindBin(layer,chInLayer), parameters[3]);
2270
2271     hnoiseTriggers->SetBinContent(hnoiseTriggers->FindBin(layer,chInLayer), ithSpectraTrigg->second.GetHG()->GetEntries());
2272   }
2273   if ( debug > 0)
2274     std::cout << "============================== done fitting" << std::endl;
2275
2276   RootOutput->cd();
2277
2278   if (IsCalibSaveToFile()){
2279     TString fileCalibPrint = RootOutputName;
2280     fileCalibPrint         = fileCalibPrint.ReplaceAll(".root","_calib.txt");
2281     calib.PrintCalibToFile(fileCalibPrint);
2282   }
2283
2284
2285   TcalibOut->Fill();
2286   TcalibOut->Write();
2287
2288   RootOutput->Write();
2289   RootOutput->Close();
2290
2291
2292   RootOutputHist->cd("IndividualCellsTrigg");
2293     for(ithSpectra=hSpectraTrigg.begin(); ithSpectra!=hSpectraTrigg.end(); ++ithSpectra){
2294       ithSpectra->second.Write(true);
2295     }
2296   RootOutputHist->cd();
2297
2298     hMeanPedHGvsCellID->Write();
2299     hMeanPedLGvsCellID->Write();
2300     hspectraHGMeanVsLayer->Write();
2301     hspectraHGSigmaVsLayer->Write();
2302     hspectraLGMeanVsLayer->Write();
2303     hspectraLGSigmaVsLayer->Write();
2304     hspectraHGvsCellID->Write();
2305     hspectraLGvsCellID->Write();
2306
2307     hnoiseTriggers->Write();
2308   // close open root files
2309   RootOutputHist->Write();
2310   RootOutputHist->Close();
2311
2312   RootInput->Close();
2313
2314   // Get run info object
2315   std::map<int,RunInfo>::iterator it=ri.find(runNr);
2316
2317   // create directory for plot output
2318   TString outputDirPlots = GetPlotOutputDir();
2319   gSystem->Exec("mkdir -p "+outputDirPlots);
2320
2321   //**********************************************************************
2322   // Create canvases for channel overview plotting
2323   //**********************************************************************
2324   Double_t textSizeRel = 0.035;
2325   StyleSettingsBasics("pdf");
2326   SetPlotStyle();
2327
2328   TCanvas* canvas2DCorr = new TCanvas("canvasCorrPlots","",0,0,1450,1300);  // gives the page size
2329   DefaultCancasSettings( canvas2DCorr, 0.08, 0.13, 0.045, 0.07);
2330
2331   canvas2DCorr->SetLogz(0);
2332   PlotSimple2D( canvas2DCorr, hspectraHGMeanVsLayer, -10000, -10000, textSizeRel, Form("%s/HG_NoiseMean.pdf", outputDirPlots.Data()), it->second, 1, kFALSE, "colz", true );
2333   PlotSimple2D( canvas2DCorr, hspectraHGSigmaVsLayer, -10000, -10000, textSizeRel, Form("%s/HG_NoiseSigma.pdf", outputDirPlots.Data()), it->second, 1, kFALSE, "colz", true);
2334   PlotSimple2D( canvas2DCorr, hspectraLGMeanVsLayer, -10000, -10000, textSizeRel, Form("%s/LG_NoiseMean.pdf", outputDirPlots.Data()), it->second, 1, kFALSE, "colz", true);
2335   PlotSimple2D( canvas2DCorr, hspectraLGSigmaVsLayer, -10000, -10000, textSizeRel, Form("%s/LG_NoiseSigma.pdf", outputDirPlots.Data()), it->second, 1, kFALSE, "colz", true);
2336   canvas2DCorr->SetLogz(1);
2337   PlotSimple2D( canvas2DCorr, hspectraHGvsCellID, -10000, -10000, textSizeRel, Form("%s/HG_Noise.pdf", outputDirPlots.Data()), it->second, 1, kFALSE, "colz", true);
2338   PlotSimple2D( canvas2DCorr, hspectraLGvsCellID, -10000, -10000, textSizeRel, Form("%s/LG_Noise.pdf", outputDirPlots.Data()), it->second, 1, kFALSE, "colz", true);
2339
2340   PlotSimple2D( canvas2DCorr, hnoiseTriggers, -10000, -10000, textSizeRel, Form("%s/LG_Noise.pdf", outputDirPlots.Data()), it->second, 1, kFALSE, "colz", true, Form( "evt. coll = %d", evts));
2341   //***********************************************************************************************************
2342   //********************************* 8 Panel overview plot  **************************************************
2343   //***********************************************************************************************************
2344   //*****************************************************************
2345     // Test beam geometry (beam coming from viewer)
2346     //===========================================================
2347     //||    8 (4)    ||    7 (5)   ||    6 (6)   ||    5 (7)   ||  row 0
2348     //===========================================================
2349     //||    1 (0)    ||    2 (1)   ||    3 (2)   ||    4 (3)   ||  row 1
2350     //===========================================================
2351     //    col 0     col 1       col 2     col  3
2352     // rebuild pad geom in similar way (numbering -1)
2353   //*****************************************************************
2354   TCanvas* canvas8Panel;
2355   TPad* pad8Panel[8];
2356   Double_t topRCornerX[8];
2357   Double_t topRCornerY[8];
2358   Int_t textSizePixel = 30;
2359   Double_t relSize8P[8];
2360   CreateCanvasAndPadsFor8PannelTBPlot(canvas8Panel, pad8Panel,  topRCornerX, topRCornerY, relSize8P, textSizePixel);
2361
2362   for (Int_t l = 0; l < setup->GetNMaxLayer()+1; l++){
2363     PlotNoiseAdvWithFitsFullLayer (canvas8Panel,pad8Panel, topRCornerX, topRCornerY, relSize8P, textSizePixel,
2364                                     hSpectra, hSpectraTrigg, setup, true, 0, 450, 1.2, l, 0,
2365                                     Form("%s/NoiseTrigg_HG_Layer%02d.pdf" ,outputDirPlots.Data(), l), it->second);
2366   }
2367
2368
2369   return true;
2370 }
2371
2372
2373
2374 //***********************************************************************************************
2375 //*********************** Calibration routine ***************************************************
2376 //***********************************************************************************************
2377 bool Analyses::Calibrate(void){
2378   std::cout<<"Calibrate"<<std::endl;
2379
2380   std::map<int,RunInfo> ri=readRunInfosFromFile(RunListInputName.Data(),debug,0);
2381
2382   std::cout << "Additional Output with histos being created: " << RootOutputNameHist.Data() << std::endl;
2383   if(Overwrite){
2384     std::cout << "recreating file with hists" << std::endl;
2385     RootOutputHist = new TFile(RootOutputNameHist.Data(),"RECREATE");
2386   } else{
2387     std::cout << "newly creating file with hists" << std::endl;
2388     RootOutputHist = new TFile(RootOutputNameHist.Data(),"CREATE");
2389   }
2390
2391   // create HG and LG histo's per channel
2392   TH2D* hspectraHGCorrvsCellID      = new TH2D( "hspectraHGCorr_vsCellID","ADC spectrum High Gain corrected vs CellID; cell ID; ADC_{HG} (arb. units)  ; counts ",
2393                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5, 4000,-200,3800);
2394   hspectraHGCorrvsCellID->SetDirectory(0);
2395   TH2D* hspectraLGCorrvsCellID      = new TH2D( "hspectraLGCorr_vsCellID","ADC spectrum Low Gain corrected vs CellID; cell ID; ADC_{LG} (arb. units)  ; counts  ",
2396                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5, 4000,-200,3800);
2397   hspectraLGCorrvsCellID->SetDirectory(0);
2398   TH2D* hspectraHGvsCellID      = new TH2D( "hspectraHG_vsCellID","ADC spectrum High Gain vs CellID; cell ID; ADC_{HG} (arb. units)   ; counts ",
2399                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5, 4000,0,4000);
2400   hspectraHGvsCellID->SetDirectory(0);
2401   TH2D* hspectraLGvsCellID      = new TH2D( "hspectraLG_vsCellID","ADC spectrum Low Gain vs CellID; cell ID; ADC_{LG} (arb. units)  ; counts",
2402                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5, 4000,0,4000);
2403   hspectraLGvsCellID->SetDirectory(0);
2404   TH2D* hspectraEnergyvsCellID  = new TH2D( "hspectraEnergy_vsCellID","Energy vs CellID; cell ID; E (mip eq./tile)    ; counts",
2405                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5, 6000,0,200);
2406   hspectraEnergyvsCellID->SetDirectory(0);
2407   TH2D* hspectraEnergyTotvsNCells  = new TH2D( "hspectraTotEnergy_vsNCells","Energy vs CellID; N_{Cells}; E_{tot} (mip eq./tile) ; counts",
2408                                             setup->GetMaxCellID()+1, -0.5, setup->GetMaxCellID()+1-0.5, 6000,0,1000);
2409   hspectraEnergyTotvsNCells->SetDirectory(0);
2410
2411   Int_t runNr = -1;
2412   RootOutput->cd();
2413   std::cout << "starting to run calibration: " << TcalibIn <<  "\t" << TcalibIn->GetEntry(0) << std::endl;
2414   TcalibIn->GetEntry(0);
2415   int actCh1st = 0;
2416   double averageScale = calib.GetAverageScaleHigh(actCh1st);
2417   double avLGHGCorr   = calib.GetAverageLGHGCorr();
2418   std::cout << "average HG mip: " << averageScale << "\t active ch: "<< actCh1st<< std::endl;
2419
2420   // setup local trigger sel
2421   TRandom3* rand    = new TRandom3();
2422   Int_t localTriggerTiles = 4;
2423   if (yearData == 2023){
2424     localTriggerTiles = 6;
2425   }
2426   double factorMinTrigg   = 0.8;
2427   double factorMaxTrigg   = 2.;
2428   if (yearData == 2023){
2429     factorMinTrigg    = 0.9;
2430     factorMaxTrigg    = 2.;
2431   }
2432
2433
2434   double minMipFrac = 0.3;
2435   int corrHGADCSwap = 3500;
2436   int evts=TdataIn->GetEntries();
2437   for(int i=0; i<evts; i++){
2438     TdataIn->GetEntry(i);
2439     if (i%5000 == 0 && debug > 0) std::cout << "Reading " <<  i << " / " << evts << " events" << std::endl;
2440     if (i == 0)runNr = event.GetRunNumber();
2441     double Etot = 0;
2442     int nCells  = 0;
2443     for(int j=0; j<event.GetNTiles(); j++){
2444       Caen* aTile=(Caen*)event.GetTile(j);
2445       // remove bad channels from output
2446       if (calib.GetBCCalib() && calib.GetBadChannel(aTile->GetCellID())!= 3 ){
2447         event.RemoveTile(aTile);
2448         j--;
2449         continue;
2450       }
2451
2452       double energy = 0;
2453       double corrHG = aTile->GetADCHigh()-calib.GetPedestalMeanH(aTile->GetCellID());
2454       double corrLG = aTile->GetADCLow()-calib.GetPedestalMeanL(aTile->GetCellID());
2455       if(corrHG<corrHGADCSwap){
2456         if(corrHG/calib.GetScaleHigh(aTile->GetCellID()) > minMipFrac){
2457           energy=corrHG/calib.GetScaleHigh(aTile->GetCellID());
2458         }
2459       } else {
2460         energy=corrLG/calib.GetCalcScaleLow(aTile->GetCellID());
2461       }
2462       hspectraHGvsCellID->Fill(aTile->GetCellID(), aTile->GetADCHigh());
2463       hspectraLGvsCellID->Fill(aTile->GetCellID(), aTile->GetADCLow());
2464       hspectraHGCorrvsCellID->Fill(aTile->GetCellID(), corrHG);
2465       hspectraLGCorrvsCellID->Fill(aTile->GetCellID(), corrLG);
2466       if(energy!=0){
2467         // calculate trigger primitives
2468         aTile->SetLocalTriggerPrimitive(event.CalculateLocalMuonTrigg(calib, rand, aTile->GetCellID(), localTriggerTiles, avLGHGCorr));
2469         bool localMuonTrigg = event.InspectIfLocalMuonTrigg(aTile->GetCellID(), averageScale, factorMinTrigg, factorMaxTrigg);
2470
2471         aTile->SetE(energy);
2472         aTile->SetLocalTriggerBit(0);
2473
2474         if (localMuonTrigg) aTile->SetLocalTriggerBit(1);
2475         hspectraEnergyvsCellID->Fill(aTile->GetCellID(), energy);
2476         Etot=Etot+energy;
2477         nCells++;
2478       } else {
2479         event.RemoveTile(aTile);
2480         j--;
2481       }
2482     }
2483     hspectraEnergyTotvsNCells->Fill(nCells,Etot);
2484     RootOutput->cd();
2485     TdataOut->Fill();
2486   }
2487   TdataOut->Write();
2488   TsetupIn->CloneTree()->Write();
2489
2490   if (IsCalibSaveToFile()){
2491     TString fileCalibPrint = RootOutputName;
2492     fileCalibPrint         = fileCalibPrint.ReplaceAll(".root","_calib.txt");
2493     calib.PrintCalibToFile(fileCalibPrint);
2494   }
2495
2496   TcalibOut->Fill();
2497   TcalibOut->Write();
2498
2499
2500   RootOutput->Close();
2501   RootInput->Close();
2502
2503   RootOutputHist->cd();
2504
2505   hspectraHGvsCellID->Write();
2506   hspectraLGvsCellID->Write();
2507   hspectraHGCorrvsCellID->Write();
2508   hspectraLGCorrvsCellID->Write();
2509   hspectraEnergyvsCellID->Write();
2510   hspectraEnergyTotvsNCells->Write();
2511
2512   TH1D* hspectraEnergyTot = (TH1D*)hspectraEnergyTotvsNCells->ProjectionY();
2513   hspectraEnergyTot->SetDirectory(0);
2514   TH1D* hspectraNCells = (TH1D*)hspectraEnergyTotvsNCells->ProjectionX();
2515   hspectraNCells->SetDirectory(0);
2516   hspectraEnergyTot->Write("hTotEnergy");
2517   hspectraNCells->Write("hNCells");
2518
2519   RootOutputHist->Close();
2520   //**********************************************************************
2521   //********************* Plotting ***************************************
2522   //**********************************************************************
2523   // Get run info object
2524   std::map<int,RunInfo>::iterator it=ri.find(runNr);
2525
2526   TString outputDirPlots = GetPlotOutputDir();
2527   gSystem->Exec("mkdir -p "+outputDirPlots);
2528
2529   //**********************************************************************
2530   // Create canvases for channel overview plotting
2531   //**********************************************************************
2532   Double_t textSizeRel = 0.035;
2533   StyleSettingsBasics("pdf");
2534   SetPlotStyle();
2535
2536   TCanvas* canvas2DCorr = new TCanvas("canvasCorrPlots","",0,0,1450,1300);  // gives the page size
2537   DefaultCancasSettings( canvas2DCorr, 0.08, 0.13, 0.045, 0.07);
2538   canvas2DCorr->SetLogz(1);
2539   PlotSimple2D( canvas2DCorr, hspectraHGvsCellID, -10000, -10000, textSizeRel, Form("%s/HG.pdf", outputDirPlots.Data()), it->second, 1, kFALSE, "colz", true);
2540   PlotSimple2D( canvas2DCorr, hspectraLGvsCellID, -10000, -10000, textSizeRel, Form("%s/LG.pdf", outputDirPlots.Data()), it->second, 1, kFALSE, "colz", true);
2541   PlotSimple2D( canvas2DCorr, hspectraHGCorrvsCellID, -10000, -10000, textSizeRel, Form("%s/HGCorr.pdf", outputDirPlots.Data()), it->second, 1, kFALSE, "colz", true);
2542   PlotSimple2D( canvas2DCorr, hspectraLGCorrvsCellID, -10000, -10000, textSizeRel, Form("%s/LGCorr.pdf", outputDirPlots.Data()), it->second, 1, kFALSE, "colz", true);
2543   PlotSimple2D( canvas2DCorr, hspectraEnergyvsCellID, -10000, -10000, textSizeRel, Form("%s/EnergyVsCellID.pdf", outputDirPlots.Data()), it->second, 1, kFALSE, "colz", true);
2544   PlotSimple2D( canvas2DCorr, hspectraEnergyTotvsNCells, -10000, -10000, textSizeRel, Form("%s/EnergyTotalVsNCells.pdf", outputDirPlots.Data()), it->second, 1, kFALSE, "colz", true);
2545
2546   TCanvas* canvas1DSimple = new TCanvas("canvas1DSimple","",0,0,1450,1300);  // gives the page size
2547   DefaultCancasSettings( canvas1DSimple, 0.08, 0.03, 0.03, 0.07);
2548   hspectraEnergyTot->Scale(1./evts);
2549   hspectraEnergyTot->GetYaxis()->SetTitle("counts/event");
2550   PlotSimple1D(canvas1DSimple, hspectraEnergyTot, -10000, -10000, textSizeRel, Form("%s/EnergyTot.pdf", outputDirPlots.Data()), it->second, 1, Form("#LT E_{Tot} #GT = %.1f (mip/tile eq.)",hspectraEnergyTot->GetMean() ));
2551   hspectraNCells->Scale(1./evts);
2552   hspectraNCells->GetYaxis()->SetTitle("counts/event");
2553   PlotSimple1D(canvas1DSimple, hspectraNCells, -10000, -10000, textSizeRel, Form("%s/NCells.pdf", outputDirPlots.Data()), it->second, 1, Form("#LT N_{Cells} #GT = %.1f",hspectraNCells->GetMean() ));
2554
2555   return true;
2556 }
2557
2558
2559 //***********************************************************************************************
2560 //*********************** Save Noise triggers only ***************************************************
2561 //***********************************************************************************************
2562 bool Analyses::SaveNoiseTriggersOnly(void){
2563   std::cout<<"Save noise triggers into separate file"<<std::endl;
2564   TcalibIn->GetEntry(0);
2565   int evts=TdataIn->GetEntries();
2566   for(int i=0; i<evts; i++){
2567     TdataIn->GetEntry(i);
2568     if (i%5000 == 0 && i > 0 && debug > 0) std::cout << "Reading " <<  i << " / " << evts << " events" << std::endl;
2569     for(int j=0; j<event.GetNTiles(); j++){
2570       Caen* aTile=(Caen*)event.GetTile(j);
2571       // testing for noise trigger
2572       if(aTile->GetLocalTriggerBit()!= (char)2){
2573         event.RemoveTile(aTile);
2574         j--;
2575       }
2576     }
2577     RootOutput->cd();
2578     TdataOut->Fill();
2579   }
2580   TdataOut->Write();
2581   TsetupIn->CloneTree()->Write();
2582
2583   if (IsCalibSaveToFile()){
2584     TString fileCalibPrint = RootOutputName;
2585     fileCalibPrint         = fileCalibPrint.ReplaceAll(".root","_calib.txt");
2586     calib.PrintCalibToFile(fileCalibPrint);
2587   }
2588
2589   TcalibOut->Fill();
2590   TcalibOut->Write();
2591   RootOutput->Close();
2592   RootInput->Close();
2593
2594   return true;
2595 }
2596
2597 //***********************************************************************************************
2598 //*********************** Save local muon triggers only ***************************************************
2599 //***********************************************************************************************
2600 bool Analyses::SaveMuonTriggersOnly(void){
2601   std::cout<<"Save local muon triggers into separate file"<<std::endl;
2602   TcalibIn->GetEntry(0);
2603   int evts=TdataIn->GetEntries();
2604   for(int i=0; i<evts; i++){
2605     TdataIn->GetEntry(i);
2606     if (i%5000 == 0 && i > 0 && debug > 0) std::cout << "Reading " <<  i << " / " << evts << " events" << std::endl;
2607     for(int j=0; j<event.GetNTiles(); j++){
2608       Caen* aTile=(Caen*)event.GetTile(j);
2609       // testing for muon trigger
2610       if(aTile->GetLocalTriggerBit()!= (char)1){
2611         event.RemoveTile(aTile);
2612         j--;
2613       }
2614     }
2615     RootOutput->cd();
2616     TdataOut->Fill();
2617   }
2618   TdataOut->Write();
2619   TsetupIn->CloneTree()->Write();
2620
2621   if (IsCalibSaveToFile()){
2622     TString fileCalibPrint = RootOutputName;
2623     fileCalibPrint         = fileCalibPrint.ReplaceAll(".root","_calib.txt");
2624     calib.PrintCalibToFile(fileCalibPrint);
2625   }
2626
2627   TcalibOut->Fill();
2628   TcalibOut->Write();
2629   RootOutput->Close();
2630   RootInput->Close();
2631
2632   return true;
2633 }
2634
2635
2636
2637 //***********************************************************************************************
2638 //*********************** Create output files ***************************************************
2639 //***********************************************************************************************
2640 bool Analyses::CreateOutputRootFile(void){
2641   if(Overwrite){
2642     RootOutput=new TFile(RootOutputName.Data(),"RECREATE");
2643   } else{
2644     RootOutput = new TFile(RootOutputName.Data(),"CREATE");
2645   }
2646   if(RootOutput->IsZombie()){
2647     std::cout<<"Error opening '"<<RootOutput<<"'no reachable path? Exist without force mode to overwrite?..."<<std::endl;
2648     return false;
2649   }
2650   return true;
2651 }
2652
2653 //***********************************************************************************************
2654 //*********************** Read external bad channel map *****************************************
2655 //***********************************************************************************************
2656 std::map<int,short> Analyses::ReadExternalBadChannelMap(void){
2657
2658   std::cout << "Reading in external mapping file" << std::endl;
2659   std::map<int,short> bcmap;
2660
2661   std::ifstream bcmapFile;
2662   bcmapFile.open(ExternalBadChannelMap,std::ios_base::in);
2663   if (!bcmapFile) {
2664     std::cout << "ERROR: file " << ExternalBadChannelMap.Data() << " not found!" << std::endl;
2665     return bcmap;
2666   }
2667
2668   for( TString tempLine; tempLine.ReadLine(bcmapFile, kTRUE); ) {
2669     // check if line should be considered
2670     if (tempLine.BeginsWith("%") || tempLine.BeginsWith("#")){
2671       continue;
2672     }
2673     if (debug > 1) std::cout << tempLine.Data() << std::endl;
2674
2675     // Separate the string according to tabulators
2676     TObjArray *tempArr  = tempLine.Tokenize(" ");
2677     if(tempArr->GetEntries()<2){
2678       if (debug > 1) std::cout << "nothing to be done" << std::endl;
2679       delete tempArr;
2680       continue;
2681     }
2682
2683     int mod     = ((TString)((TObjString*)tempArr->At(0))->GetString()).Atoi();
2684     int layer   = ((TString)((TObjString*)tempArr->At(1))->GetString()).Atoi();
2685     int row     = ((TString)((TObjString*)tempArr->At(2))->GetString()).Atoi();
2686     int col     = ((TString)((TObjString*)tempArr->At(3))->GetString()).Atoi();
2687     short bc    = short(((TString)((TObjString*)tempArr->At(4))->GetString()).Atoi());
2688
2689     int cellID  = setup->GetCellID( row, col, layer, mod);
2690
2691     if (debug > 1) std::cout << "cellID " << cellID << "\t BC status: " << bc<< std::endl;
2692     bcmap[cellID]=bc;
2693   }
2694   std::cout << "registered " << bcmap.size() << " bad channels!" << std::endl;
2695   return bcmap;
2696
2697 }