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