Back to home page

sPhenix code displayed by LXR

 
 

    


File indexing completed on 2025-08-05 08:20:27

0001 // use #include "" only for your local include and put
0002 // those in the first line(s) before any #include <>
0003 // otherwise you are asking for weird behavior
0004 // (more info - check the difference in include path search when using "" versus <>)
0005 
0006 #include "LocalPolMon.h"
0007 
0008 #include <onlmon/OnlMon.h>  // for OnlMon
0009 #include <onlmon/OnlMonDB.h>
0010 #include <onlmon/OnlMonServer.h>
0011 
0012 #include <Event/msg_profile.h>
0013 #include <caloreco/CaloWaveformFitting.h>
0014 
0015 #include <Event/Event.h>
0016 #include <Event/EventTypes.h>
0017 #include <Event/eventReceiverClient.h>
0018 #include <Event/msg_profile.h>
0019 
0020 #include <TH1.h>
0021 #include <TH2.h>
0022 #include <TObjArray.h>
0023 #include <TObjString.h>
0024 #include <TProfile.h>
0025 #include <TRandom.h>
0026 #include <TString.h>
0027 #include <TSystem.h>
0028 
0029 #include <cmath>
0030 #include <cstdio>  // for printf
0031 #include <ctime>
0032 #include <fstream>
0033 #include <iostream>
0034 #include <map>
0035 #include <sstream>
0036 #include <string>  // for allocator, string, char_traits
0037 #include <utility>
0038 
0039 LocalPolMon::LocalPolMon(const std::string& name)
0040   : OnlMon(name)
0041 {
0042   // leave ctor fairly empty, its hard to debug if code crashes already
0043   // during a new LocalPolMon()
0044   return;  // test
0045 }
0046 
0047 LocalPolMon::~LocalPolMon()
0048 {
0049   // you can delete NULL pointers it results in a NOOP (No Operation)
0050   return;
0051 }
0052 
0053 int LocalPolMon::Init()
0054 {
0055   iPoint = 0;
0056   for (bool& i : goodtrigger)
0057   {  // 16 triggers for gl1p
0058     i = false;
0059   }
0060   goodtrigger[10] = true;
0061   OnlMonServer* se = OnlMonServer::instance();
0062   const char* locpolcal = getenv("LOCALPOLCALIB");
0063   if (!locpolcal)
0064   {
0065     std::cout << "LOCALPOLCALIB environment variable not set" << std::endl;
0066     exit(1);
0067   }
0068   std::string lpconfigfilename = std::string(locpolcal) + "/" + "LocPolConfig.dat";
0069   std::ifstream config(lpconfigfilename);
0070   std::ostringstream msg_config(lpconfigfilename);
0071 
0072   if (!config)
0073   {
0074     msg_config << lpconfigfilename << " could not be opened.";
0075     se->send_message(this, MSG_SOURCE_LOCALPOL, MSG_SEV_FATAL, msg_config.str(), 2);
0076     exit(1);
0077   }
0078   int err_counter = 0;
0079   TString key;
0080   TString val;
0081   while (!config.eof())
0082   {
0083     config >> key >> val;
0084     std::cout << "Config in " << config.good() << " state\n";
0085     if (!config.good())
0086     {
0087       break;
0088     }
0089     key.ToLower();
0090     std::cout << "Reading " << key << " with value " << val << std::endl;
0091     if (key == "trigger")
0092     {
0093       if (!val.IsDigit())
0094       {
0095         std::cout << key << ": expecting an integer (2^(trigger bit number))\n keep only bit 10 as default value.\n";
0096         goodtrigger[10] = true;
0097       }
0098       else
0099       {
0100         int ival = val.Atoi();
0101         if (ival < 0 || ival > 65535)
0102         {
0103           std::cout << key << ": value outside expected range [0; 65535]\n Keep only bit 10 as default value.\n";
0104           goodtrigger[10] = true;
0105         }
0106         else
0107         {
0108           for (int i = 0; i < 16; i++)
0109           {
0110             goodtrigger[i] = false;
0111             if ((ival & (1 << i)))
0112             {
0113               goodtrigger[i] = true;
0114             }
0115           }
0116         }
0117       }
0118     }
0119     else if (key == "testfake")
0120     {
0121       if (!val.IsDigit())
0122       {
0123         std::cout << key << ": expecting 0/1 for true or false\n Keep false as default.\n";
0124         fake = false;
0125       }
0126       else
0127       {
0128         int ival = val.Atoi();
0129         if (ival != 0 && ival != 1)
0130         {
0131           std::cout << key << ": expecting 0/1 for true or false\n Keep false as default.\n";
0132           fake = false;
0133         }
0134         else
0135         {
0136           if (ival == 1)
0137           {
0138             fake = true;
0139             erc = nullptr;
0140           }
0141           else
0142           {
0143             fake = false;
0144           }
0145         }
0146       }
0147     }
0148     else if (key == "monitoring"){
0149       val.ToLower();
0150       if (val == "online"){
0151     erc = new eventReceiverClient("gl1daq");
0152       }
0153       else if (val == "offline" ){
0154     erc = new eventReceiverClient("localhost");
0155       }
0156       else {
0157     std::cout<< key << ": expecting either online (data stream)/offline (with eventServer -d 5263250 -s 5 -i -v -f path to gl1 prdf) "<<std::endl;
0158     std::cout<<"Fall back to online monitoring"<<std::endl;
0159     erc = new eventReceiverClient("gl1daq");
0160       }
0161     }
0162     else if (key == "sphenixgap")
0163     {
0164       if (!val.IsDigit())
0165       {
0166         std::cout << key << ": expecting an integer as the first bunch number of the continuous gap sequence\n Keep 111 as default.\n";
0167         ExpectedsPhenixGapPosition = 111;
0168       }
0169       else
0170       {
0171         int ival = val.Atoi();
0172         if (ival >= 0 && ival < 120)
0173         {
0174           ExpectedsPhenixGapPosition = ival;
0175         }
0176         else
0177         {
0178           std::cout << key << ": value outside expected range [0;119]\n Keep 117 as default value\n";
0179           ExpectedsPhenixGapPosition = 111;
0180         }
0181       }
0182     }
0183     else if (key == "threshold")
0184     {
0185       if (!val.IsDigit())
0186       {
0187         std::cout << key << ": expecting a positive integer\n Keep 6000 as default value\n";
0188         EventCountThresholdGap = 6000;
0189       }
0190       else
0191       {
0192         int ival = val.Atoi();
0193         if (ival > 0 && ival < 1e5)
0194         {
0195           EventCountThresholdGap = ival;
0196         }
0197         else
0198         {
0199           std::cout << key << ": integer seems too big (>1e5)\n Keep default 6000\n";
0200           EventCountThresholdGap = 6000;
0201         }
0202       }
0203     }
0204     else if (key == "verbosity")
0205     {
0206       if (!val.IsDigit())
0207       {
0208         std::cout << key << ": expecting 0/1 for false or true\n Keep 0(false) as default\n";
0209         verbosity = false;
0210       }
0211       else
0212       {
0213         int ival = val.Atoi();
0214         if (ival != 0 && ival != 1)
0215         {
0216           std::cout << key << ": value should be 0 or 1 for false of true\n Keep 0(false) as default\n";
0217           verbosity = false;
0218         }
0219         else if (ival == 1)
0220         {
0221           std::cout << "Making it verbose" << std::endl;
0222           verbosity = true;
0223         }
0224       }
0225     }
0226     else if (key == "integrated")
0227     {
0228       if (!val.IsDigit())
0229       {
0230         std::cout << key << ": expecting 0/1 for false or true\n Keep 0(false) as default\n";
0231         Integrated = false;
0232       }
0233       else
0234       {
0235         int ival = val.Atoi();
0236         if (ival != 0 && ival != 1)
0237         {
0238           std::cout << key << ": value should be 0 or 1 for false of true\n Keep 0(false) as default\n";
0239           Integrated = false;
0240         }
0241         else if (ival == 1)
0242         {
0243           std::cout << "Making it integrated" << std::endl;
0244           Integrated = true;
0245         }
0246       }
0247     }
0248     else if (key == "x0north")
0249     {
0250       if (val.IsFloat())
0251       {
0252         if (val.Atof() < -3. || val.Atof() > 3.)
0253         {
0254           std::cout << key << ": value outside the expected range [-3.;3]cm\n keeping 0.0 as default value." << std::endl;
0255         }
0256         else
0257         {
0258           ZeroPosition[1] = val.Atof();
0259         }
0260       }
0261       else
0262       {
0263         std::cout << key << ": value is not a float\n keeping 0.0 as default value" << std::endl;
0264         ZeroPosition[1] = 0.;
0265       }
0266     }
0267     else if (key == "y0north")
0268     {
0269       if (val.IsFloat())
0270       {
0271         if (val.Atof() < -3. || val.Atof() > 3.)
0272         {
0273           std::cout << key << ": value outside the expected range [-3.;3]cm\n keeping 0.0 as default value." << std::endl;
0274         }
0275         else
0276         {
0277           ZeroPosition[0] = val.Atof();
0278         }
0279       }
0280       else
0281       {
0282         std::cout << key << ": value is not a float\n keeping 0.0 as default value" << std::endl;
0283         ZeroPosition[0] = 0.;
0284       }
0285     }
0286     else if (key == "x0south")
0287     {
0288       if (val.IsFloat())
0289       {
0290         if (val.Atof() < -3. || val.Atof() > 3.)
0291         {
0292           std::cout << key << ": value outside the expected range [-3.;3]cm\n keeping 0.0 as default value." << std::endl;
0293         }
0294         else
0295         {
0296           ZeroPosition[3] = val.Atof();
0297         }
0298       }
0299       else
0300       {
0301         std::cout << key << ": value is not a float\n keeping 0.0 as default value" << std::endl;
0302         ZeroPosition[3] = 0.;
0303       }
0304     }
0305     else if (key == "y0south")
0306     {
0307       if (val.IsFloat())
0308       {
0309         if (val.Atof() < -3. || val.Atof() > 3.)
0310         {
0311           std::cout << key << ": value outside the expected range [-3.;3]cm\n keeping 0.0 as default value." << std::endl;
0312         }
0313         else
0314         {
0315           ZeroPosition[2] = val.Atof();
0316         }
0317       }
0318       else
0319       {
0320         std::cout << key << ": value is not a float\n keeping 0.0 as default value" << std::endl;
0321         ZeroPosition[2] = 0.;
0322       }
0323     }
0324     else if (key == "thresholdasymnewpoint")
0325     {
0326       if (!val.IsDigit())
0327       {
0328         std::cout << key << ": value is expected to be a positive integer\n Keeping default value of 5000 events." << std::endl;
0329       }
0330       else
0331       {
0332         std::cout << key << "Changed from " << EventsAsymmetryNewPoint << " to " << val.Atoi() << std::endl;
0333         EventsAsymmetryNewPoint = val.Atoi();
0334       }
0335     }
0336     else if (key == "zdc1cut")
0337       {
0338     if(!val.IsFloat())
0339       {
0340         std::cout<< key << " : value is expected to be a positive float\n Keeping default value of 60."<<std::endl; 
0341       }
0342     else
0343       {
0344         std::cout<< key << " Changed from "<< ZDC1Cut << " to "<<val.Atof()<<std::endl;
0345         ZDC1Cut = val.Atof();
0346       }
0347       }
0348     else if (key == "zdc2cut")
0349       {
0350     if(!val.IsFloat())
0351       {
0352         std::cout<< key << " : value is expected to be a positive float\n Keeping default value of 15."<<std::endl; 
0353       }
0354     else
0355       {
0356         std::cout<< key << " Changed from "<< ZDC2Cut << " to "<<val.Atof()<<std::endl;
0357         ZDC2Cut = val.Atof();
0358       }
0359       }
0360     else if (key == "vetocut")
0361       {
0362     if(!val.IsFloat())
0363       {
0364         std::cout<< key << " : value is expected to be a positive float\n Keeping default value of 150."<<std::endl; 
0365       }
0366     else
0367       {
0368         std::cout<< key << " Changed from "<< VetoCut << " to "<<val.Atof()<<std::endl;
0369         VetoCut = val.Atof();
0370       }
0371       }
0372     else if (key == "multiplicitylow")
0373       {
0374     if(!val.IsDigit())
0375       {
0376         std::cout<< key << " : value is expected to be a positive integer\n Keeping default value of 1"<<std::endl; 
0377       }
0378     else
0379       {
0380         std::cout<< key << " Changed from "<< MultiLow << " to "<<val.Atoi()<<std::endl;
0381         MultiLow = val.Atoi();
0382       }
0383       }
0384     else if (key == "multiplicityhigh")
0385       {
0386     if(!val.IsDigit())
0387       {
0388         std::cout<< key << " : value is expected to be a positive integer\n Keeping default value of 7 (for vertical and +1 for horizontal)"<<std::endl; 
0389       }
0390     else
0391       {
0392         std::cout<< key << " Changed from "<< MultiHigh << " to "<<val.Atoi()<<std::endl;
0393         MultiHigh = val.Atoi();
0394       }
0395       }
0396     else if (key == "smdthreshold")
0397       {
0398     if(!val.IsFloat())
0399       {
0400         std::cout<< key << " : value is expected to be a positive float\n Keeping default value of 5"<<std::endl; 
0401       }
0402     else
0403       {
0404         std::cout<< key << " Changed from "<< SMDthr << " to "<<val.Atof()<<std::endl;
0405         SMDthr = val.Atof();
0406       }
0407       }
0408     else
0409     {
0410       err_counter++;
0411       std::cout << "Unknown configuration \n Expected:\n -verbosity\n -threshold\n -sphenixgap\n -testfake\n -monitoring\n -multiplicitylow\n -multiplicityHigh\n -trigger\n -X0north\n -X0south\n -Y0north\n -Y0south\n -thresholdasymnewpoint\n -nVetocut\n -SMDthreshold\n -ZDC1cut\n -ZDC2cut\n -Integrated\n key words" << std::endl;
0412     }
0413     if (err_counter > 3)
0414     {
0415       std::cout << "More than 3 unknown key words, abort" << std::endl;
0416       exit(1);
0417     }
0418   }
0419 
0420   std::string lpmappingfilename = std::string(locpolcal) + "/" + "ChannelMapping.txt";
0421   std::ifstream mapping(lpmappingfilename);
0422   std::ostringstream msg_mapping(lpmappingfilename);
0423 
0424   if (!mapping)
0425   {
0426     msg_mapping << lpmappingfilename << " could not be opened.";
0427     se->send_message(this, MSG_SOURCE_LOCALPOL, MSG_SEV_FATAL, msg_mapping.str(), 2);
0428     exit(1);
0429   }
0430   int adc, array, lowcut, highcut;
0431   std::string ChannelName;
0432   for (int i = 0; i < 128; i++)
0433   {
0434     mapping >> adc >> array >> lowcut >> highcut >> ChannelName;
0435     Chmapping[adc] = array;
0436     if(array>-1){
0437       lowSample[array]=lowcut;
0438       highSample[array]=highcut;
0439     }
0440   }
0441 
0442   const char* zdccalib = getenv("ZDCCALIB");
0443   if (!zdccalib)
0444   {
0445     std::cout << "ZDCCALIB environment variable not set" << std::endl;
0446     exit(1);
0447   }
0448   std::string calibfilename = std::string(zdccalib) + "/" + "SMDRelativeGain.dat";
0449   std::ifstream calibinput(calibfilename);
0450   std::ostringstream msg(calibfilename);
0451 
0452   if (!calibinput)
0453   {
0454     msg << calibfilename << " could not be opened.";
0455     se->send_message(this, MSG_SOURCE_LOCALPOL, MSG_SEV_FATAL, msg.str(), 2);
0456     exit(1);
0457   }
0458   int NS;
0459   int channel;
0460   float gain;
0461   for (int ch = 0; ch < 32; ch++)
0462   {
0463     calibinput >> NS >> channel >> gain;
0464 
0465     if (NS == 0 && channel >= 0 && channel < 16)
0466     {
0467       smd_north_relatgain[channel] = gain;
0468     }
0469     else if (NS == 1 && channel >= 0 && channel < 16)
0470     {
0471       smd_south_relatgain[channel] = gain;
0472     }
0473     else
0474     {
0475       msg << calibfilename << "could not understand the content, expect int(0/1) int(0-15) float for North or South, channel number and relative gain";
0476       se->send_message(this, MSG_SOURCE_ZDC, MSG_SEV_FATAL, msg.str(), 2);
0477       exit(1);
0478     }
0479   }
0480 
0481   h_Counts = new TH1D*[4];
0482   h_CountsScramble = new TH1D*[4];
0483 
0484   h_Counts[0] = new TH1D("h_BlueCountsUD", "h_BlueCountsUD", 4, 0, 4);
0485   h_Counts[1] = new TH1D("h_BlueCountsLR", "h_BlueCountsLR", 4, 0, 4);
0486   h_Counts[2] = new TH1D("h_YellCountsUD", "h_YellCountsUD", 4, 0, 4);
0487   h_Counts[3] = new TH1D("h_YellCountsLR", "h_YellCountsLR", 4, 0, 4);
0488   h_CountsScramble[0] = new TH1D("h_BlueCountsScrambleUD", "h_BlueCountsScrambleUD", 4, 0, 4);
0489   h_CountsScramble[1] = new TH1D("h_BlueCountsScrambleLR", "h_BlueCountsScrambleLR", 4, 0, 4);
0490   h_CountsScramble[2] = new TH1D("h_YellCountsScrambleUD", "h_YellCountsScrambleUD", 4, 0, 4);
0491   h_CountsScramble[3] = new TH1D("h_YellCountsScrambleLR", "h_YellCountsScrambleLR", 4, 0, 4);
0492 
0493   for(int i=0; i<4; i++){
0494     se->registerHisto(this,h_Counts[i]);
0495     se->registerHisto(this,h_CountsScramble[i]);
0496   }
0497 
0498   h_time = new TProfile("h_times", "h_time", 5000, -0.5, 4999.5);
0499   se->registerHisto(this, h_time);
0500 
0501   TString BeamName[2] = {"Blue", "Yell"};
0502   TString MethodName[2] = {"Arithmetic", "Geometric"};
0503   //TString Orientation[2] = {"LR", "UD"};
0504   TString Orientation[2] = {"UD", "LR"};
0505   h_Asym = new TH1D***[2];
0506   h_AsymScramble = new TH1D***[2];
0507   for (int beam = 0; beam < 2; beam++)
0508   {
0509     h_Asym[beam] = new TH1D**[2];
0510     h_AsymScramble[beam] = new TH1D**[2];
0511     for (int method = 0; method < 2; method++)
0512     {
0513       h_Asym[beam][method] = new TH1D*[2];
0514       h_AsymScramble[beam][method] = new TH1D*[2];
0515       for (int orient = 0; orient < 2; orient++)
0516       {
0517         h_Asym[beam][method][orient] = new TH1D(Form("h_Asym%s%s%s", BeamName[beam].Data(), MethodName[method].Data(), Orientation[orient].Data()), Form("Fwd %s %s %s Asym.", BeamName[beam].Data(), Orientation[orient].Data(), MethodName[method].Data()), 5000, -0.5, 4999.5);
0518         h_AsymScramble[beam][method][orient] = new TH1D(Form("h_AsymScramble%s%s%s", BeamName[beam].Data(), MethodName[method].Data(), Orientation[orient].Data()), Form("Bwk %s %s %s Asym.", BeamName[beam].Data(), Orientation[orient].Data(), MethodName[method].Data()), 5000, -0.5, 4999.5);
0519 
0520         se->registerHisto(this, h_Asym[beam][method][orient]);
0521         se->registerHisto(this, h_AsymScramble[beam][method][orient]);
0522       }
0523     }
0524   }
0525   for(int itrig=0; itrig<16; itrig++){
0526     h_trigger[itrig]=new TH1D(Form("h_trigger%d",itrig),Form("h_trigger%d",itrig),120,0,120);
0527     se->registerHisto(this,h_trigger[itrig]);
0528   }
0529   h_events=new TH1D("h_events","h_events",20,0,20);
0530   se->registerHisto(this,h_events);
0531 
0532   hspinpattern=new TH2I("hspinpattern","hspinpattern",120,0,120,2,0,2);
0533   se->registerHisto(this,hspinpattern);
0534 
0535   hmultiplicity[0]=new TH1D("hmultiplicitySMD_NH","hmultiplicitySMD_NH",9,0,9);
0536   hmultiplicity[1]=new TH1D("hmultiplicitySMD_NV","hmultiplicitySMD_NV",8,0,8);
0537   hmultiplicity[2]=new TH1D("hmultiplicitySMD_SH","hmultiplicitySMD_SH",9,0,9);
0538   hmultiplicity[3]=new TH1D("hmultiplicitySMD_SV","hmultiplicitySMD_SV",8,0,8);
0539 
0540   hposition[0]=new TH1D("hpositionSMD_NH_up","hpositionSMD_NH_up",100,-5,5);
0541   hposition[1]=new TH1D("hpositionSMD_NV_up","hpositionSMD_NV_up",100,-5,5);
0542   hposition[2]=new TH1D("hpositionSMD_SH_up","hpositionSMD_SH_up",100,-5,5);
0543   hposition[3]=new TH1D("hpositionSMD_SV_up","hpositionSMD_SV_up",100,-5,5);
0544   hposition[4]=new TH1D("hpositionSMD_NH_dn","hpositionSMD_NH_dn",100,-5,5);
0545   hposition[5]=new TH1D("hpositionSMD_NV_dn","hpositionSMD_NV_dn",100,-5,5);
0546   hposition[6]=new TH1D("hpositionSMD_SH_dn","hpositionSMD_SH_dn",100,-5,5);
0547   hposition[7]=new TH1D("hpositionSMD_SV_dn","hpositionSMD_SV_dn",100,-5,5);
0548 
0549   hadcsum[0]=new TH1D("hadcsumSMD_NH","hadcsumSMD_NH",100,0,4.5);
0550   hadcsum[1]=new TH1D("hadcsumSMD_NV","hadcsumSMD_NV",100,0,4.5);
0551   hadcsum[2]=new TH1D("hadcsumSMD_SH","hadcsumSMD_SH",100,0,4.5);
0552   hadcsum[3]=new TH1D("hadcsumSMD_SV","hadcsumSMD_SV",100,0,4.5);
0553 
0554   for(int i=0; i<4; i++){
0555     se->registerHisto(this,hmultiplicity[i]);
0556     se->registerHisto(this,hposition[i]);
0557     se->registerHisto(this,hposition[i+4]);
0558     se->registerHisto(this,hadcsum[i]);
0559   }
0560 
0561   hwaveform[0]=new TH2D("hwaveform0","hwaveform0",16,-0.5,15.5,16384,0,16384);
0562   hwaveform[1]=new TH2D("hwaveform1","hwaveform1",16,-0.5,15.5,16384,0,16384);
0563   hwaveform[2]=new TH2D("hwaveform2","hwaveform2",16,-0.5,15.5,16384,0,16384);
0564   hwaveform[3]=new TH2D("hwaveform3","hwaveform3",16,-0.5,15.5,16384,0,16384);
0565   hwaveform[4]=new TH2D("hwaveform4","hwaveform4",16,-0.5,15.5,16384,0,16384);
0566   hwaveform[5]=new TH2D("hwaveform5","hwaveform5",16,-0.5,15.5,16384,0,16384);
0567   for(int i=0; i<6; i++){
0568     se->registerHisto(this,hwaveform[i]);
0569   }
0570   
0571   //htimesync=new TH2I("htimesync","htimesync",65000,0,65000,65000,0,65000);
0572   //se->registerHisto(this,htimesync);
0573   hsyncfrac=new TH1D("hsyncfrac","hsyncfrac",2,0,2);
0574   se->registerHisto(this,hsyncfrac);
0575 
0576   Bluespace=new TH2D("Bluespace","Bluespace",50,-5,5,50,-5,5);
0577   se->registerHisto(this,Bluespace);
0578  
0579   Yellowspace=new TH2D("Yellowspace","Yellowspace",50,-5,5,50,-5,5);
0580   se->registerHisto(this,Yellowspace);
0581 
0582   hclocks = new TH2D("hclocks","hclocks",8192,0,8192,8192,0,8192);
0583   se->registerHisto(this, hclocks);
0584 
0585   hevolsync = new TH2D("hevolsync","",10000,0,30000000,2,0,2);
0586   se->registerHisto(this,hevolsync);
0587 
0588   hshiftevol=new TH2D("hshiftevol","",10000,0,30000000,20,0,20);
0589   se->registerHisto(this,hshiftevol);
0590   
0591   WaveformProcessingFast = new CaloWaveformFitting();
0592   myRandomBunch = new TRandom(0);
0593   Reset();
0594   return 0;
0595 }
0596 
0597 int LocalPolMon::BeginRun(const int runno)
0598 {
0599   // if you need to read calibrations on a run by run basis
0600   // this is the place to do it
0601 
0602   // Initialisation of the map (hopefully this step will not be required when gl1p scalers become available
0603   for (int i = 0; i < 16; i++)
0604   {  // 16 triggers for gl1p
0605     for (int bunchnr = 0; bunchnr < 120; bunchnr++)
0606     {
0607       gl1_counter[i][bunchnr] = 0;
0608     }
0609     if (verbosity)
0610     {
0611       std::cout << goodtrigger[15 - i];
0612     }
0613   }
0614   if (verbosity)
0615   {
0616     std::cout << "\n";
0617   }
0618   EvtShift=0;
0619   Initfirstbunch=false;
0620   Reset();
0621   RetrieveSpinPattern(runno);
0622   // goodtrigger[10]=true;//In 2023, it was MBD N&S
0623   return 0;
0624 }
0625 
0626 float LocalPolMon::anaWaveformFast(Packet* p, const int channel, const int low, const int high, const int ihisto)
0627 {
0628   std::vector<float> waveform;
0629   waveform.reserve(high-low);
0630   //waveform.reserve(p->iValue(0, "SAMPLES"));
0631   //for (int s = 0; s < p->iValue(0, "SAMPLES"); s++)
0632   for (int s = low; s < high; s++)
0633   {
0634     waveform.push_back(p->iValue(s, channel));
0635     hwaveform[ihisto]->Fill(s,p->iValue(s, channel));
0636   }
0637   std::vector<std::vector<float>> multiple_wfs;
0638   multiple_wfs.push_back(waveform);
0639 
0640   std::vector<std::vector<float>> fitresults_zdc;
0641   fitresults_zdc = WaveformProcessingFast->calo_processing_fast(multiple_wfs);
0642 
0643   std::vector<float> result;
0644   result = fitresults_zdc.at(0);
0645   return result.at(0);
0646 }
0647 
0648 int LocalPolMon::process_event(Event* e /* evt */)
0649 {
0650 
0651   //if (e->getEvtType() == BEGRUNEVENT /*9*/)
0652   //{  // spin patterns stored in BeginRun event
0653   //  RetrieveSpinPattern(e);
0654   //}
0655 
0656   RetrieveTriggerDistribution(e);
0657   h_events->Fill(0);
0658   /******Here we do the matching to identify the crossing shift   ******************/
0659   if (evtcnt > EventCountThresholdGap){
0660     StartAbortGapData = RetrieveAbortGapData();
0661   }
0662   else
0663   {
0664     StartAbortGapData = ExpectedsPhenixGapPosition;  // default value from config
0665   }
0666 
0667   CrossingShift = StartAbortGapPattern - StartAbortGapData;
0668   if (verbosity)
0669   {
0670     std::cout << "Crossing shift: " << CrossingShift << " = ( " << StartAbortGapPattern << " - " << StartAbortGapData << " )\n";
0671   }
0672   /************* end of crossing shift *********************/
0673 
0674   Packet* psmd = e->getPacket(packetid_smd);
0675   int bunchnr = 0;
0676   if (psmd)
0677   {
0678     h_events->Fill(1);
0679     long long int zdc_clock=psmd->lValue(0,"CLOCK");
0680     bunchnr=RetrieveBunchNumber(e,zdc_clock);
0681     //Did not manage to retrieve the proper bunch number, no need to continue the calculation
0682     //let us clean and leave
0683     if(bunchnr<0){
0684       delete psmd;
0685       psmd=nullptr;
0686       return 0;
0687     }
0688     h_events->Fill(2);
0689 
0690     if (fake)
0691     {
0692       bunchnr += 15 + myRandomBunch->Integer(4);
0693     }
0694 
0695     // get minimum on ZDC second module
0696     // get minimum on ZDC second module
0697     signalZDCN1 = anaWaveformFast(psmd, ZDCN1, lowSample[ZDCN1], highSample[ZDCN1],0);
0698     signalZDCS1 = anaWaveformFast(psmd, ZDCS1, lowSample[ZDCS1], highSample[ZDCS1],1);
0699     signalZDCN2 = anaWaveformFast(psmd, ZDCN2, lowSample[ZDCN2], highSample[ZDCN2],0);
0700     signalZDCS2 = anaWaveformFast(psmd, ZDCS2, lowSample[ZDCS2], highSample[ZDCS2],1);
0701     vetoNF = anaWaveformFast(psmd, 16, lowSample[ivetoNF], highSample[ivetoNF],4);
0702     vetoNB = anaWaveformFast(psmd, 17, lowSample[ivetoNB], highSample[ivetoNB],4);
0703     vetoSF = anaWaveformFast(psmd, 80, lowSample[ivetoSF], highSample[ivetoSF],5);
0704     vetoSB = anaWaveformFast(psmd, 81, lowSample[ivetoSB], highSample[ivetoSB],5);
0705 
0706     if ( (signalZDCN2 < 10 || signalZDCN1<75) && (signalZDCS2 < 10 || signalZDCS1<75) )
0707     {
0708       delete psmd;
0709       psmd=nullptr;
0710       return 0;
0711     }
0712     h_events->Fill(3);
0713 
0714     // for(int ch=16; ch<47; ch++){//according to mapping in ZDC logbook entry #85 March 9th 2024
0715     for (auto& it : Chmapping)
0716     {  // new mapping for ZDC/SMD/Veto with 2 ADC boards//May 13th 2024
0717       if (it.second < 16)
0718       {
0719         continue;
0720       }
0721       if (it.second > 47)
0722       {
0723         continue;
0724       }
0725       float signalFast = anaWaveformFast(psmd, it.first, lowSample[it.second],highSample[it.second],it.second/32+2);  // fast waveform fitting
0726 
0727       int ch = it.second;
0728       // Scale according to relative gain calibration factor
0729       //if (ch < 16){
0730       //    zdc_adc[ch]=signalFast;
0731       //}
0732       //else 
0733       if (ch < 32)
0734       {
0735         smd_adc[ch - 16] = signalFast * smd_north_relatgain[ch - 16];
0736       }
0737       else
0738       {
0739         smd_adc[ch - 16] = signalFast * smd_south_relatgain[ch - 32];
0740       }
0741     }
0742     float Weights[4] = {0};
0743     memset(Weights, 0, sizeof(Weights));
0744     float AveragePosition[4] = {0};
0745     memset(AveragePosition, 0, sizeof(AveragePosition));
0746 
0747     for (int ch = 0; ch < 8; ch++)
0748     {
0749       Weights[0] += smd_adc[ch];
0750       AveragePosition[0] += ConversionSign[0] * pitchY * (ch - (nchannelsY-1.) / 2) * smd_adc[ch];  // North Y direction (Asym in blue)
0751       Weights[2] += smd_adc[ch + 16];
0752       AveragePosition[2] += ConversionSign[2] * pitchY * (ch - (nchannelsY-1.) / 2) * smd_adc[ch + 16];  // South Y direction (Asym in Yellow)
0753 
0754       if (ch == 7)
0755       {
0756         continue;
0757       }
0758       Weights[1] += smd_adc[ch + 8];
0759       AveragePosition[1] += ConversionSign[1] * pitchX * (ch - (nchannelsX-1.) / 2) * smd_adc[ch + 8];  // North X direction (Asym in Blue)
0760       Weights[3] += smd_adc[ch + 24];
0761       AveragePosition[3] += ConversionSign[3] * pitchX * (ch - (nchannelsX-1.) / 2) * smd_adc[ch + 24];  // South X direction (Asym in Yellow)
0762     }
0763 
0764     for (int i = 0; i < 4; i++)
0765     {
0766       if (Weights[i] > 0.0)
0767       {
0768         AveragePosition[i] /= Weights[i];
0769     h_events->Fill(4+i);
0770       }
0771       else
0772       {
0773     AveragePosition[i]=0.;
0774         continue;  // most likely the most appropriate
0775       }
0776 
0777       if (AveragePosition[i] < ZeroPosition[i] - 0.5)
0778       {
0779     h_events->Fill(8+i);
0780     if (verbosity)
0781       {
0782         std::cout<<"Spin pattern size for "<<i/2<<" : "<<SpinPatterns[i/2].size()<<std::endl;
0783       }
0784         // (i/2)=0 for blue beam=North, (i/2)=1 for yellow beam=South
0785         if (SpinPatterns[i / 2].at((120 + bunchnr + CrossingShift) % 120) > 0)
0786         {
0787           if(GoodSelection(i/2)) {
0788         h_Counts[i]->Fill(3);  // Right for pointing up
0789         hposition[i]->Fill(AveragePosition[i]);
0790       }
0791         }
0792         else if (SpinPatterns[i / 2].at((120 + bunchnr + CrossingShift) % 120) < 0)
0793         {
0794           if(GoodSelection(i/2)) {
0795         h_Counts[i]->Fill(1);  // Left for pointing down
0796         hposition[i+4]->Fill(AveragePosition[i]);
0797       }
0798         }
0799 
0800         // we swap the spin pattern to get random orientation and check biased asymmetry
0801         if (SpinPatterns[i / 2 == 0 ? 1 : 0].at((120 + bunchnr + CrossingShift) % 120) > 0)
0802         {
0803           if(GoodSelection(i/2))h_CountsScramble[i]->Fill(3);  // Right for pointing up
0804         }
0805         else if (SpinPatterns[i / 2 == 0 ? 1 : 0].at((120 + bunchnr + CrossingShift) % 120) < 0)
0806         {
0807           if(GoodSelection(i/2))h_CountsScramble[i]->Fill(1);  // Left for pointing down
0808         }
0809       }
0810       else if (AveragePosition[i] > ZeroPosition[i] + 0.5)
0811       {
0812     h_events->Fill(8+i);
0813         // (i/2)=0 for blue beam, =1 for yellow beam
0814         if (SpinPatterns[i / 2].at((120 + bunchnr + CrossingShift) % 120) > 0)
0815         {
0816           if(GoodSelection(i/2)){
0817         h_Counts[i]->Fill(0);  // Left for pointing up
0818         hposition[i]->Fill(AveragePosition[i]);
0819       }
0820         }
0821         else if (SpinPatterns[i / 2].at((120 + bunchnr + CrossingShift) % 120) < 0)
0822         {
0823           if(GoodSelection(i/2)){
0824         h_Counts[i]->Fill(2);  // Right for pointing down
0825         hposition[i+4]->Fill(AveragePosition[i]);
0826       }
0827         }
0828 
0829         // we swap the spin pattern to get random orientation and check biased asymmetry
0830         if (SpinPatterns[i / 2 == 0 ? 1 : 0].at((120 + bunchnr + CrossingShift) % 120) > 0)
0831         {
0832           if(GoodSelection(i/2))h_CountsScramble[i]->Fill(0);  // Left for pointing up
0833         }
0834         else if (SpinPatterns[i / 2 == 0 ? 1 : 0].at((120 + bunchnr + CrossingShift) % 120) < 0)
0835         {
0836           if(GoodSelection(i/2))h_CountsScramble[i]->Fill(2);  // Right for pointing down
0837         }
0838       }
0839     }
0840     double rnorth2=pow(AveragePosition[0]-ZeroPosition[0],2)+pow(AveragePosition[1]-ZeroPosition[1],2);
0841     double rsouth2=pow(AveragePosition[2]-ZeroPosition[2],2)+pow(AveragePosition[3]-ZeroPosition[3],2);
0842     if(rnorth2>1 && GoodSelection(0)){
0843       Bluespace->Fill(AveragePosition[1],AveragePosition[0]);
0844     }
0845     if(rsouth2>1 && GoodSelection(1)){
0846       Yellowspace->Fill(AveragePosition[3],AveragePosition[2]);
0847     }
0848     h_time->Fill(iPoint, e->getTime());
0849     evtcnt++;
0850     evtcntA++;
0851     delete psmd;
0852     psmd=nullptr;
0853   }
0854   else
0855   {
0856     if (verbosity)
0857     {
0858       std::cout << "Missing ZDC/SMD packet " << packetid_smd << std::endl;
0859     }
0860   }
0861 
0862   /**** Compute asymmetries if we have enough events ****************/
0863   //if (evtcnt % EventsAsymmetryNewPoint == 0)
0864   if (evtcntA > EventsAsymmetryNewPoint )
0865   {
0866     h_events->Fill(12);
0867     for (int ibeam = 0; ibeam < 2; ibeam++)
0868     {
0869       for (int orient = 0; orient < 2; orient++)
0870       {
0871         double L_U = h_Counts[2 * ibeam + orient]->GetBinContent(1);
0872         double R_D = h_Counts[2 * ibeam + orient]->GetBinContent(2);
0873         double L_D = h_Counts[2 * ibeam + orient]->GetBinContent(3);
0874         double R_U = h_Counts[2 * ibeam + orient]->GetBinContent(4);
0875         // if(L_U<0){
0876         //   std::cout<<iPoint<<" "<<evtcnt<<" "<<h_Counts[2*ibeam+orient]->GetEntries()<<" "<<std::endl;
0877         // }
0878         double* asymresult = ComputeAsymmetries(L_U, R_D, L_D, R_U);
0879 
0880         h_Asym[ibeam][0][orient]->SetBinContent(iPoint + 1, asymresult[0]);
0881         h_Asym[ibeam][0][orient]->SetBinError(iPoint + 1, asymresult[1]);
0882 
0883         h_Asym[ibeam][1][orient]->SetBinContent(iPoint + 1, asymresult[2]);
0884         h_Asym[ibeam][1][orient]->SetBinError(iPoint + 1, asymresult[3]);
0885 
0886         L_U = h_CountsScramble[2 * ibeam + orient]->GetBinContent(1);
0887         R_D = h_CountsScramble[2 * ibeam + orient]->GetBinContent(2);
0888         L_D = h_CountsScramble[2 * ibeam + orient]->GetBinContent(3);
0889         R_U = h_CountsScramble[2 * ibeam + orient]->GetBinContent(4);
0890 
0891         double* scrambleresult = ComputeAsymmetries(L_U, R_D, L_D, R_U);
0892 
0893         h_AsymScramble[ibeam][0][orient]->SetBinContent(iPoint + 1, scrambleresult[0]);
0894         h_AsymScramble[ibeam][0][orient]->SetBinError(iPoint + 1, scrambleresult[1]);
0895 
0896         h_AsymScramble[ibeam][1][orient]->SetBinContent(iPoint + 1, scrambleresult[2]);
0897         h_AsymScramble[ibeam][1][orient]->SetBinError(iPoint + 1, scrambleresult[3]);
0898 
0899         delete asymresult;
0900         delete scrambleresult;
0901       }
0902     }
0903     // Now that everything has been calculated, let move to the next point for next event
0904     iPoint++;
0905     evtcntA=0;
0906     if(!Integrated){
0907       for(int i =0; i<4; i++){
0908     h_Counts[i]->Reset();
0909     h_CountsScramble[i]->Reset();
0910       }
0911     }
0912   }
0913 
0914   return 0;
0915 }
0916 
0917 int LocalPolMon::Reset()
0918 {
0919   // reset our internal counters
0920   evtcnt = 0;
0921   evtcntA=0;
0922   iPoint = 0;
0923   failuredepth=0;
0924   for (int i = 0; i < 4; i++)
0925   {
0926     h_Counts[i]->Reset();
0927     h_CountsScramble[i]->Reset();
0928   }
0929   h_time->Reset();
0930   return 0;
0931 }
0932 
0933 bool LocalPolMon::GoodSelection(int i)
0934 {
0935   bool goodselection=true;
0936   if(i==0){//dealing with north side
0937     int nv=0; 
0938     int nh=0;
0939     double sumH=0;
0940     double sumV=0;
0941     for(int ch=0; ch<8; ch++) {
0942       nh +=(smd_adc[ch]  >SMDthr)?1:0;
0943       sumH+=smd_adc[ch];
0944     }
0945     for(int ch=0; ch<7; ch++) {
0946       nv +=(smd_adc[ch+8]>SMDthr)?1:0;
0947       sumV+=smd_adc[ch+8];
0948     }
0949     goodselection &=(signalZDCN1>ZDC1Cut && signalZDCN2>ZDC2Cut);
0950     goodselection &=(vetoNF<VetoCut&&vetoNB<VetoCut);
0951     if(goodselection&&(nv>MultiLow&&nv<MultiHigh+1)) hmultiplicity[0]->Fill(nh);
0952     if(goodselection&&(nh>MultiLow&&nh<MultiHigh+1)) hmultiplicity[1]->Fill(nv);
0953     goodselection &=(nv>MultiLow&&nh>MultiLow);
0954     goodselection &=(nv<MultiHigh&&nh<MultiHigh+1);
0955     if((sumH>0)&&goodselection) hadcsum[0]->Fill(log10(sumH));
0956     if((sumV>0)&&goodselection) hadcsum[1]->Fill(log10(sumV));
0957   }
0958   else{//south side
0959     int nv=0; 
0960     int nh=0;
0961     double sumH=0;
0962     double sumV=0;
0963     for(int ch=0; ch<8; ch++) {
0964       nh +=(smd_adc[ch+16] >SMDthr)?1:0;
0965       sumH+=smd_adc[ch+16];
0966     }
0967     for(int ch=0; ch<7; ch++) {
0968       nv +=(smd_adc[ch+24] >SMDthr)?1:0;
0969       sumV+=smd_adc[ch+24];
0970     }
0971     goodselection &=(signalZDCS1>ZDC1Cut && signalZDCS2>ZDC2Cut);
0972     goodselection &=(vetoSF<VetoCut&&vetoSB<VetoCut);
0973     if(goodselection&&(nv>MultiLow&&nv<MultiHigh+1)) hmultiplicity[2]->Fill(nh);
0974     if(goodselection&&(nh>MultiLow&&nh<MultiHigh+1)) hmultiplicity[3]->Fill(nv);
0975     goodselection &=(nv>MultiLow&&nh>MultiLow);
0976     goodselection &=(nv<MultiHigh+1&&nh<MultiHigh+1);
0977     if((sumH>0)&&goodselection) hadcsum[2]->Fill(log10(sumH));
0978     if((sumV>0)&&goodselection) hadcsum[3]->Fill(log10(sumV));
0979   }
0980 
0981   return goodselection;
0982 }
0983 
0984 double* LocalPolMon::ComputeAsymmetries(double L_U, double R_D, double L_D, double R_U)
0985 {
0986   double* result = new double[4];
0987   double leftA = L_U + R_D;
0988   double rightA = L_D + R_U;
0989   double tmpNumA = leftA - rightA;
0990   double tmpDenA = leftA + rightA;
0991   result[0] = 0;
0992   result[1] = 0;
0993 
0994   if (tmpDenA > 0)
0995   {
0996     result[0] = tmpNumA / tmpDenA;
0997     result[1] = 2 * sqrt(pow(rightA, 2) * leftA + pow(leftA, 2) * rightA) / pow(tmpDenA, 2);
0998   }
0999   //std::cout<<leftA<<" "<<rightA<<"\t\t"<<result[0]<<" +/- "<<result[1]<<std::endl;
1000 
1001   double leftG = sqrt(L_U * R_D);
1002   double rightG = sqrt(L_D * R_U);
1003   double tmpNumG = leftG - rightG;
1004   double tmpDenG = leftG + rightG;
1005   result[2] = 0;
1006   result[3] = 0;
1007 
1008   if (tmpDenG > 0)
1009   {
1010     result[2] = tmpNumG / tmpDenG;
1011     result[3] = sqrt(pow(rightG, 2) * leftA + pow(leftG, 2) * rightA) / pow(tmpDenG, 2);
1012   }
1013   return result;
1014 }
1015 
1016 void LocalPolMon::RetrieveSpinPattern(int runnb){
1017   SpinPatterns[BLUE].clear();
1018   SpinPatterns[YELLOW].clear();
1019   for(int i = 0; i < 120; i++){
1020     SpinPatterns[BLUE][i]=0;
1021     SpinPatterns[YELLOW][i]=0;
1022   }
1023   TString retvalBlue;
1024   TString retvalYellow;
1025   const char* user=getenv("USER");
1026   if(strcmp(user,"phnxrc")==0){
1027     //We are on machine which can access operator1 computer without password
1028     retvalBlue=gSystem->GetFromPipe(Form("$ONLMON_MACROS/GetSpinPatternHTTP.sh %d",14902));
1029     retvalYellow=gSystem->GetFromPipe(Form("$ONLMON_MACROS/GetSpinPatternHTTP.sh %d",14903));
1030     if(verbosity){
1031       std::cout<<"Retrieving SpinPattern from HTTP"<<std::endl;
1032     }
1033   }
1034   else{
1035     OnlMonServer* se = OnlMonServer::instance();
1036     TString runtype=se->GetRunType();
1037     runtype.ToLower();
1038     retvalBlue=gSystem->GetFromPipe(Form("$ONLMON_MACROS/GetSpinPatternGL1.sh %d %s %d",14902,runtype.Data(),runnb));
1039     retvalYellow=gSystem->GetFromPipe(Form("$ONLMON_MACROS/GetSpinPatternGL1.sh %d %s %d",14903,runtype.Data(),runnb));
1040     if(verbosity){
1041       std::cout<<"Retrieving SpinPattern from GL1 "<<runtype.Data() <<" prdf file for run "<<runnb<<std::endl;
1042     }
1043   }
1044   int flip=-1;//https://phenix-intra.sdcc.bnl.gov/phenix/WWW/offline/wikioff/index.php/Spin_Database_:_Note_for_Analyzer
1045   TObjArray* BunchSpinBlue=retvalBlue.Tokenize(" ");
1046   BunchSpinBlue->SetOwner(kTRUE);
1047   int nFilledBunchesBlue=BunchSpinBlue->GetEntries();
1048   if(verbosity){
1049     std::cout<<"Blue spin patterns for "<<nFilledBunchesBlue<<" bunches "<<std::endl;
1050     std::cout<<retvalBlue<<std::endl;
1051   }
1052   int stepBlue=(112-nFilledBunchesBlue)/28+1;
1053   //for(int i=0; i<120; i=i+stepBlue){
1054   for(int i=0; i<nFilledBunchesBlue; i++){
1055     SpinPatterns[BLUE][stepBlue*i]=flip*((TObjString*)BunchSpinBlue->At(i))->String().Atoi();
1056     hspinpattern->Fill(stepBlue*i,1,((TObjString*)BunchSpinBlue->At(i))->String().Atoi());
1057   }
1058   BunchSpinBlue->Clear();
1059   delete BunchSpinBlue;
1060 
1061   TObjArray* BunchSpinYellow=retvalYellow.Tokenize(" ");
1062   BunchSpinYellow->SetOwner(kTRUE);
1063   int nFilledBunchesYell=BunchSpinYellow->GetEntries();
1064   if(verbosity){
1065     std::cout<<"Yellow spin patterns for "<<nFilledBunchesYell<<" bunches "<<std::endl;
1066     std::cout<<retvalYellow<<std::endl;
1067   }
1068   int stepYellow=(112-nFilledBunchesYell)/28+1;
1069   //for(int i=0; i<120; i=i+stepYellow){
1070   for(int i=0; i<nFilledBunchesYell; i++){
1071     SpinPatterns[YELLOW][stepYellow*i]=flip*((TObjString*)BunchSpinYellow->At(i))->String().Atoi();
1072     hspinpattern->Fill(stepYellow*i,0.,((TObjString*)BunchSpinYellow->At(i))->String().Atoi());
1073   }
1074   //Check if consitent information (as the method works only with multiples of 28 filled bunches and a constant abort gap of 9 bunches
1075   if(stepYellow==stepBlue){
1076     nEmptyAbort=stepYellow-1;
1077   }
1078   BunchSpinYellow->Clear();
1079   delete BunchSpinYellow;
1080 }
1081 
1082 void LocalPolMon::RetrieveTriggerDistribution(Event* e){
1083   Event* egl1p=nullptr;
1084   Packet* pgl1p = nullptr;
1085   if (erc){
1086     if(verbosity){
1087       std::cout<<"Inside RetrieveTrigger::ERC"<<e->getEvtSequence()<<" "<<EvtShift<<std::endl;
1088     }
1089     egl1p = erc->getEvent(e->getEvtSequence()+EvtShift);
1090     if(verbosity){
1091       std::cout<<egl1p<<std::endl;
1092     }
1093   }
1094   if (egl1p){
1095     pgl1p = egl1p->getPacket(packetid_gl1);
1096     if (pgl1p){
1097       int bunchnr = pgl1p->lValue(0, "BunchNumber");
1098       for (int i=0; i<16; i++){//auto& i : gl1_counter)
1099     // 16 triggers for gl1p
1100     // With prdf pgl1->lValue(i,2); simply returns the current processed event number (which can shaddow the abort gap: the lagging bunch# at some point get back to the position of the others)
1101     // So instead, we increment the number of processed events per bunch number, for the various triggers
1102     //gl1_counter[i][bunchnr]+= (pgl1p->lValue(i,2)>0)?1:0;
1103     if(pgl1p->lValue(i,2)>0){
1104       gl1_counter[i][bunchnr]++;
1105       h_trigger[i]->Fill(bunchnr);
1106     }
1107       }
1108       delete pgl1p;
1109       pgl1p=nullptr;
1110     }
1111     delete egl1p;
1112     egl1p=nullptr;
1113   }
1114 }
1115 
1116 int LocalPolMon::RetrieveAbortGapData(){
1117   // Here we do the magic to identify the abort gap
1118   std::vector<int> gap[16];
1119   std::map<int, int> begingap;
1120   std::map<int, int> endgap;
1121 
1122   for (int i = 0; i < 16; i++){
1123     if (!goodtrigger[i]){
1124       continue;
1125     }
1126     if(h_trigger[i]->GetEntries()<1000) continue;
1127     std::map<int, long long> tmpmap = gl1_counter[i];
1128     for (int emptyfill = 0; emptyfill < 9; emptyfill++){
1129       int myminimum = min_element(tmpmap.begin(), tmpmap.end(), [](const std::pair<int, long long>& lhs, const std::pair<int, long long>& rhs)
1130                   { return lhs.second < rhs.second; })
1131     ->first;
1132       if (myminimum < 111){
1133     gap[i].push_back(120 + myminimum);
1134       }
1135       else{
1136     gap[i].push_back(myminimum);
1137       }
1138       tmpmap.erase(myminimum);
1139     }
1140     begingap[i] = (*min_element(gap[i].begin(), gap[i].end()))+nEmptyAbort;
1141     endgap[i] = (*max_element(gap[i].begin(), gap[i].end()));
1142     if (endgap[i] - begingap[i] > 9){
1143       std::cout << " Weird abort gap is larger than 9 bunches " << endgap[i] << " - " << begingap[i] << " for trigger " << i << std::endl;
1144     }
1145   }
1146   for (auto ib = begingap.begin(); ib != begingap.end(); ++ib){
1147     if (begingap.begin()->second != ib->second)
1148       {
1149     std::cout << " Weird abort gap not in the same location between trigger bit 0 and trigger bit " << ib->second << std::endl;
1150       }
1151   }
1152   if (begingap.empty()) return 111;
1153   else return (begingap.begin()->second) % 120;
1154 }
1155 
1156 
1157 int LocalPolMon::RetrieveBunchNumber(Event* e, long long int zdc_clock){
1158   Event* egl1 = nullptr;
1159   int bunch=-1;
1160   //static int localfail=0;
1161   if(!Initfirstbunch){
1162     EvtShift=0;
1163     EvtShiftValid=0;
1164     Prevzdc_clock=zdc_clock;
1165   }
1166   if(verbosity){
1167     std::cout<<"Inside RetrieveBunchNumber"<<std::endl;
1168   }
1169   if(failuredepth>1000){
1170     if(verbosity){
1171       std::cout<<"ZDC and GL1p asynchronous by more than 1000 events"<<std::endl;
1172       std::cout<<"Giving up this event and go to the next one"<<std::endl;
1173     }
1174     EvtShift=EvtShiftValid;
1175     failuredepth=0;
1176     return -1;
1177   }
1178   //std::cout<<e->getEvtSequence()<<std::endl;
1179   if (erc){
1180     if(verbosity){
1181       std::cout<<"Inside RetrieveBunchNumber::ERC "<<e->getEvtSequence() <<" "<<EvtShift <<std::endl;
1182     }
1183     egl1 = erc->getEvent(e->getEvtSequence()+EvtShift);
1184     if(verbosity){
1185       std::cout<<egl1<<std::endl;
1186     }
1187   }
1188   if (egl1){
1189     if(verbosity){
1190       std::cout<<"Inside RetrieveBunchNumber::GL1"<<std::endl;
1191     }
1192     Packet* pgl1p = egl1->getPacket(packetid_gl1);
1193     if (pgl1p){
1194       if(verbosity){
1195     std::cout<<"Inside RetrieveBunchNumber::PGL1p"<<std::endl;
1196       }
1197       long long int gl1_clock=pgl1p->lValue(0, "BCO");
1198       if(verbosity){
1199        std::cout<<"EvtShift: "<<EvtShift<<" zdc: "<<zdc_clock<<"    gl1p: "<<gl1_clock<<std::endl;
1200       }
1201 
1202       if(!Initfirstbunch){
1203     Prevgl1_clock=gl1_clock;
1204     Initfirstbunch=true;
1205     bunch = pgl1p->lValue(0, "BunchNumber");
1206     if(verbosity){
1207       std::cout<<"Init Bunch number is : "<<bunch<<std::endl;
1208     }
1209     delete pgl1p;
1210     pgl1p=nullptr;
1211     delete egl1;
1212     egl1=nullptr;
1213     return bunch;
1214       }
1215       if(zdc_clock<Prevzdc_clock){
1216     //Prevzdc_clock=Prevzdc_clock-4294967296;//despite the long long, it seems it is only 32 bits
1217     //std::cerr<<"Mismatched: "<< e->getEvtSequence()<<" shift: "<<EvtShift<<" zdc: "<<zdc_clock<<" - "<<Prevzdc_clock<<" = "<<(zdc_clock-Prevzdc_clock) <<"    gl1p: "<<gl1_clock<<" - "<<Prevgl1_clock<<" = "<<(gl1_clock-Prevgl1_clock) <<std::endl;
1218     zdc_clock+=(long long int)1<<32;
1219     //std::cerr<<"And now Mismatched: "<< e->getEvtSequence()<<" shift: "<<EvtShift<<" zdc: "<<zdc_clock<<" - "<<Prevzdc_clock<<" = "<<(zdc_clock-Prevzdc_clock) <<"    gl1p: "<<gl1_clock<<" - "<<Prevgl1_clock<<" = "<<(gl1_clock-Prevgl1_clock) <<std::endl;
1220     //exit(1);
1221       }
1222       hclocks->Fill((gl1_clock-Prevgl1_clock)%8192,(zdc_clock-Prevzdc_clock)%8192);
1223       if((gl1_clock-Prevgl1_clock)!=(zdc_clock-Prevzdc_clock)){
1224     if(verbosity){
1225       std::cout<<"Mismatched: "<<EvtShift<<" zdc: "<<(zdc_clock-Prevzdc_clock) <<"    gl1p: "<<(gl1_clock-Prevgl1_clock) <<std::endl;
1226         }
1227     //std::cerr<<"Mismatched: "<< e->getEvtSequence()<<" shift: "<<EvtShift<<" zdc: "<<zdc_clock<<" - "<<Prevzdc_clock<<" = "<<(zdc_clock-Prevzdc_clock) <<"    gl1p: "<<gl1_clock<<" - "<<Prevgl1_clock<<" = "<<(gl1_clock-Prevgl1_clock) <<std::endl;
1228     //if((zdc_clock-Prevzdc_clock)<0){
1229     //  exit(1) ;
1230     //}
1231     EvtShift++;
1232     delete pgl1p;
1233     pgl1p=nullptr;
1234     delete egl1;
1235     egl1=nullptr;
1236     failuredepth++;
1237     hsyncfrac->Fill(0.);
1238     hevolsync->Fill(e->getEvtSequence(),0);
1239     bunch=RetrieveBunchNumber(e,zdc_clock);
1240       }
1241       else{
1242     Prevgl1_clock=gl1_clock;
1243     if(zdc_clock>(long long int)1<<32){
1244       zdc_clock-=(long long int)1<<32;
1245     }
1246     Prevzdc_clock=zdc_clock;
1247     EvtShiftValid=EvtShiftValid;
1248     hsyncfrac->Fill(1.);
1249     hevolsync->Fill(e->getEvtSequence(),1);
1250     hshiftevol->Fill(e->getEvtSequence(),EvtShift);
1251     bunch = pgl1p->lValue(0, "BunchNumber");
1252     //failuredepth=0;
1253     if(verbosity){
1254       std::cout<<"Bunch number is : "<<bunch<<std::endl;
1255     }
1256       }
1257       delete pgl1p;
1258       pgl1p = nullptr;
1259     }
1260     else{
1261       if (verbosity){
1262     std::cout << "Failed grabing gl1 from event receiver, Bunch number unknown" << std::endl;
1263       }
1264     }
1265     delete egl1;
1266     egl1 = nullptr;
1267   }
1268   return bunch;
1269 }