onlinemonitor/rocmonitor/TFiberHodProc.cxx (r4864/r3723)

Go to the documentation of this file.

#include "TFiberHodProc.h"

#include "Riostream.h"

#include "TH1.h"
#include "TH2.h"
#include "TROOT.h"
#include "TSystem.h"
#include "TPedestalExtractor.h"

#include "TGo4MbsEvent.h"
#include "TGo4WinCond.h"
#include "TGo4Log.h"

#include "roc/Message.h"

#include <algorithm>

#include "roc/Message.h"
//#include "roc/Board.h"
#include "TRocEvent.h"
#include "TFiberHodParam.h"



TFiberHodProc::TFiberHodProc() :
   TCBMBeamtimeProc(),
   fParam(0),
   fRocEvent(0),
   fHodEvent(0),
   fHitbank(),
   fClusterbank()
{
   // dummy for tree streamer
}



//***********************************************************
// this one is used in standard factory

TFiberHodProc::TFiberHodProc(const char* name) :
   TCBMBeamtimeProc(name),
   fParam(0),
   fRocEvent(0),
   fHodEvent(0)
{
   TGo4Log::Info("TFiberHodProc: Create instance %s", name);

   TString path = GetName();
   path += "/";
   path += GetName();

   // Booking histograms
   hNrHits =       MakeTH1('I', path+"_NrHits","Hit number",        64,-0.5,63.5);
   hMaxCh =        MakeTH1('I', path+"_MaxChannel","maxChannel",    128,-0.5,127.5);
   hMaxFiberall =  MakeTH1('I', path+"_MaxFiberall","maxFiber",     128,0.5,128.5);
   hMaxFiberpl1 =  MakeTH1('I', path+"_MaxFiber1","maxFiber_pl1",   64,0.5,64.5);
   hMaxFiberpl2 =  MakeTH1('I', path+"_MaxFiber2","maxFiber_pl2",   64,0.5,64.5);
   hADC =          MakeTH1('I', path+"_ADC","ADC value",            100, 0., 4000.);
   hMaxADC =       MakeTH2('I', path+"_MaxADC","maxADC",            128,-0.5,127.5,100,0.,4000);
   hMaxADCraw =    MakeTH2('I', path+"_MaxADCraw","maxADCraw",      128,-0.5,127.5,100,0.,4000);
   hMaxFiber2D =   MakeTH2('I', path+"_Hits2D","Hits2D",            64,0.5,64.5,64,0.5,64.5);
   hNrClusters2D = MakeTH2('I', path+"_NrClusters2d","NrClusters2D",10,-0.5,9.5,10,-0.5,9.5);
   hNrClusters =   MakeTH1('I', path+"_NrClusters","NrClusters",    10,-0.5,9.5);
   hPosition2D =   MakeTH2('I', path+"_Position2D","Position2D",    64,-32*1.01,32*1.01, 64,-32*1.01,32*1.01);
   hMaxClusterTimediff = MakeTH1('I', path+"_MaxClusterTimediff", "MaxClusterTimediff", 100,-50,50);
   // Booking histograms for single event display

   fADC_cond = MakeWinCond(path+"_ADCcut", 200., 9999., hADC->GetName());

   TString parname = GetName(); parname += "Par";
   TString setupmacro = "set_"; setupmacro+=parname; setupmacro+=".C";

   // SL TODO: after next go4 release change uncomment usage of setup macro
   fParam = (TFiberHodParam *) MakeParameter(parname, "TFiberHodParam"/*, setupmacro.Data() */);

   // here we optionally override parameter values with setup macro, if existing:
   ExecuteScript(setupmacro.Data());

   // initialize and fill decoding table
   InitializeDecoding();

   fSingleEventSlice=0;
   hSingleEventHit2D = 0;
   for (Int_t icluster=1; icluster<=5; icluster++) hSingleEventCluster2D[icluster-1]=0;

   if (fParam->SingleEventMode) {
      hSingleEventHit2D = MakeTH2('F', path+"_SingleEventHit2D","SingleEventHit2D",100,0.5,100.5,128,0.5,128.5);
      for (Int_t icluster=1; icluster<=5; icluster++)
         hSingleEventCluster2D[icluster-1] =
           MakeTH2('F',Form(path+"_SingleEventCluster2D%i",icluster),
                       Form("SingleEventCluster_%i",icluster),
                       100,0.5,100.5,128,0.5,128.5);
   }

   // preserve some memory, should be enough for most applications
   fHitbank.reserve(100);
   fClusterbank.reserve(100);
}

TFiberHodProc::~TFiberHodProc()
{
   //cout << "**** TFiberHodProc: Delete instance " << endl;
}

void TFiberHodProc::InitializeDecoding(Int_t altflag)
{
   // Initialze decoding table : vector fDecoding,
   // array index is the readout channel [0..127]
   for (Int_t i=0; i<127; i++) {
      fDecoding[i].FEBchannel=i;
      fDecoding[i].fiber=-1;
      fDecoding[i].plane=-1;
      fDecoding[i].pixel=-1;
      fDecoding[i].lcn=-1;
      fDecoding[i].cable=-1;
   }

   for (Int_t ifiber=1; ifiber<=64; ifiber++) {
      // Calculate febch [0..63] for given fiber number [1..64]
      // lcn: linearconnectornumber, is the wire number on one of the
      // flat cables. [1..16]
      // each 16 fibers go to one connector.
      // fibersubnr[0..15] linear fiber counter in groups of 16

      Int_t fibersubnr=(ifiber-1)%16;

      Int_t lcn=15-fibersubnr*2;
      if (fibersubnr>=8) lcn=(fibersubnr-7)*2;

      Int_t channel=-1;
      Int_t cable=(ifiber-1)/16+1;
      Int_t pixel= ((lcn-1)/2)*8 +((lcn-1)%2);
      if (cable==1) {
         channel=(lcn-1)*4+0;
         pixel=pixel+1;
      }
      if (cable==2) {
         channel=(lcn-1)*4+2;
         pixel=pixel+3;
      }
      if (cable==3) {
         channel=(lcn-1)*4+1;
         pixel=pixel+5;
      }
      if (cable==4) {
         channel=(lcn-1)*4+3;
         pixel=pixel+7;
      }
      
      // new code to resolve cabling problem during cern-oct12
      int ifiber_bis = ifiber;
      if (fParam->WrongCabling) {
         if (ifiber <= 8 )  ifiber_bis = ifiber + 56; else
         if (ifiber <= 16 ) ifiber_bis = ifiber + 40; else
         if (ifiber <= 24 ) ifiber_bis = ifiber + 24; else
         if (ifiber <= 32 ) ifiber_bis = ifiber + 8; else
         if (ifiber <= 40 ) ifiber_bis = ifiber - 8; else
         if (ifiber <= 48 ) ifiber_bis = ifiber - 24; else
         if (ifiber <= 56 ) ifiber_bis = ifiber - 40; else
         if (ifiber <= 64 ) ifiber_bis = ifiber - 56;
         
         // and swap at the end
         ifiber_bis = 65 - ifiber_bis;
      }

      fDecoding[channel].fiber = ifiber_bis;
      fDecoding[channel].plane=1;
      fDecoding[channel].pixel=pixel;
      fDecoding[channel].lcn=lcn;
      fDecoding[channel].cable=cable;
      fDecoding[channel].position=(ifiber_bis-32.5)*1.01;

      fDecoding[channel+64].fiber = ifiber_bis;
      fDecoding[channel+64].plane=2;
      fDecoding[channel+64].pixel=pixel;
      fDecoding[channel+64].lcn=lcn;
      fDecoding[channel+64].cable=cable+4;
      fDecoding[channel+64].position=(ifiber_bis-32.5)*1.01;
   }

   TGo4Log::Info("Decoding table for Fiber hodoscope %s (rocid %i, NXid %i) initialized:",
           GetName(), fParam->RocID, fParam->NxID);
   for (Int_t i=0; i<127; i++) {
      TGo4Log::Debug("Channel %3i: FiberNr /plane: %2i /%1i, PMT PixelNr: %2i, position: %4.1f",
            fDecoding[i].FEBchannel, fDecoding[i].fiber, fDecoding[i].plane,
            fDecoding[i].pixel, fDecoding[i].position);
   }
}

void TFiberHodProc::InitEvent(TGo4EventElement* outevnt)
{
   // first assign input event:
   // since input event object is never discarded within processor lifetime,
   // we just search for subevent by name once to speed up processing
   if(fRocEvent==0) {
      TCBMBeamtimeEvent* btevent=dynamic_cast<TCBMBeamtimeEvent*>(GetInputEvent());
      if(btevent)
         fRocEvent = dynamic_cast<TRocEvent*>(btevent->GetSubEvent("ROC"));
      else
         fRocEvent = dynamic_cast<TRocEvent*>(GetInputEvent());

      if(fRocEvent==0) {
         GO4_STOP_ANALYSIS_MESSAGE("**** TFiberHodProc: Fatal error: no input event TRocEvent!!! STOP GO4");
      }
   }
   // then assign output event
   // since output event object is never discarded within processor lifetime,
   // we just search for subevent by name once to speed up processing
   if(fHodEvent==0) {
      TCBMBeamtimeEvent* btevent = dynamic_cast<TCBMBeamtimeEvent*>(outevnt);
      if(btevent)
         fHodEvent = dynamic_cast<TFiberHodEvent*>(btevent->GetSubEvent(GetName()));
      else
         fHodEvent = dynamic_cast<TFiberHodEvent*>(outevnt);

      if(fHodEvent==0) {
         GO4_STOP_ANALYSIS_MESSAGE("**** TFiberHodProc: Fatal error: output event is not a TFiberHodEvent!!! STOP GO4");
      }
   }

}

void TFiberHodProc::FinalizeEvent()
{
   if (fParam->RocID<0) return;

   if ((fRocEvent==0) || (fHodEvent==0)) {
      TGo4Log::Error("TFiberHodProc: %s ERROR: no input/output event specified", GetName());
      return;
   }

   TRocData* RocData = dynamic_cast<TRocData*>(fRocEvent->getEventElement(fParam->RocID));
   if(RocData==0) {
      TGo4Log::Error("TFiberHodProc: %s ERROR: can not find roc of id %d in input, please check configuration!!!", GetName(), fParam->RocID);
      return;
   }

   fHitbank.clear();
   fClusterbank.clear();

   //variables for finding maximum fiber and max cluster in each layer,
   Double_t maxADC[2]={-1000,-1000};
   Double_t maxADCraw[2]={-1000,-1000};
   Int_t maxCh[2]={-1,-1};
   Int_t maxFiber[2]={-1,-1};

   Int_t MaxCluster[2] = {-1, -1};
   Int_t NrClustersPerPlane[2] = { 0, 0 };  

   unsigned NrMsgs = RocData->fExtMessages.size();
   
   hNrHits->Fill(NrMsgs);

   // loop over all messages for Fiberhodoscope to fill hitbank
   // apply cut on ADCrange and TIMErange

   for (unsigned MessageNr=0; MessageNr<NrMsgs; MessageNr++) {
      TRocMessageExtended& msg = RocData->fExtMessages.at(MessageNr);

      if (msg.GetMessageType() != roc::MSG_HIT) continue;  // only interested in HIT messages
      if (msg.GetRocNumber() != fParam->RocID || msg.GetNxNumber() != fParam->NxID) continue; // only hodoscope hits

      // apply cuts on hit selection
      if (!fADC_cond->Test(msg.GetCorrectedNxADC())) continue;

      // extract message information
      hitbankentry_t newhit;

      newhit.FEBchannel = msg.GetNxChNum();   // 0..127
      newhit.fiber = fDecoding[msg.GetNxChNum()].fiber; // 1-64
      newhit.plane = fDecoding[msg.GetNxChNum()].plane; // 1,2
      newhit.ADC = msg.GetCorrectedNxADC();
      newhit.ADCraw = msg.GetNxADC();
      newhit.time = msg.GetTriggerDeltaT();
#ifdef HODOSCOPE_FULLTIMESTAMP
      newhit.fulltime = msg.GetFullTime();
#endif
      newhit.position = fDecoding[msg.GetNxChNum()].position;
      newhit.status = 0;

      fHitbank.push_back(newhit);

      // store ADC, ADCraw, maxCh, for maximum hit per each hodoscope plane.
      if (newhit.ADC > maxADC[newhit.plane-1]) {
         maxADC[newhit.plane-1] = newhit.ADC;
         maxADCraw[newhit.plane-1] = newhit.ADCraw;
         maxFiber[newhit.plane-1] = newhit.fiber;
         maxCh[newhit.plane-1] = newhit.FEBchannel;
      }
   }

   if (fParam->Debug) PrintAllHits();


   if (fHitbank.size()!=0) {

      // Do cluster finding in all hits:
      // First sort hitbank according to ADC value
      // Then look for first not-yet-clustered ADC hit and start new cluster
      // Loop over rest of hitbank, and look for neighboring hits within +-fMaxClusterDist
      // add all hits within range to the cluster, and finally calculate mean position, ADCsum
      // store cluster in clusterbank, and restart with next not-yet-clustered hit

      sort(fHitbank.begin(), fHitbank.end(), adccompare());

      // Find first not-yet-clustered hit with highest ADC
      unsigned ifirst;
      Bool_t done;
      do {
         done=kTRUE;
         for (ifirst=0; ifirst<fHitbank.size(); ifirst++) {
            if (fHitbank[ifirst].status==0) {
               done=kFALSE;
               break;
            }
         }

         if (!done) {  //Found a hit, start new cluster
            clusterbankentry_t newcluster;
            newcluster.NrHits=1;
            newcluster.plane = fHitbank[ifirst].plane;
            newcluster.meanposition = fHitbank[ifirst].ADC*fHitbank[ifirst].position;
            newcluster.meanfiber = fHitbank[ifirst].ADC*fHitbank[ifirst].fiber;
            newcluster.sumADC = fHitbank[ifirst].ADC;
            newcluster.meantime = fHitbank[ifirst].time;
#ifdef HODOSCOPE_FULLTIMESTAMP
            newcluster.meanfulltime = fHitbank[ifirst].fulltime;
#endif
            fHitbank[ifirst].status=1;      // Mark hit in hitbank as already used

            // look for additional hits for cluster, starting at position of
            // first hit in hitbank(ADC sorted)
            for (unsigned inext=ifirst+1; inext<fHitbank.size(); inext++) {
               if (fHitbank[inext].status!=0) continue;                      // hit already clustered or marked as crosstalk
               if (fHitbank[inext].plane!=fHitbank[ifirst].plane) continue;  // hit in different hodoscope plane
               // spatial distance to large
               if (abs(fHitbank[inext].fiber-fHitbank[ifirst].fiber)>fParam->MaxClusterDist) continue;
               // time difference to large
               if (fabs(fHitbank[inext].time-fHitbank[ifirst].time)>fParam->MaxTimeDiff) continue;

               // hit qualified to join the cluster
               newcluster.NrHits++;
               fHitbank[inext].status=1;
               newcluster.sumADC += fHitbank[inext].ADC;
               newcluster.meantime += fHitbank[inext].time;
#ifdef HODOSCOPE_FULLTIMESTAMP
               newcluster.meanfulltime += fHitbank[inext].fulltime;
#endif
               newcluster.meanposition += fHitbank[inext].ADC*fHitbank[inext].position;
               newcluster.meanfiber += fHitbank[inext].ADC*fHitbank[inext].fiber;
            }

            // Calculate cluster position as weighted mean, TDC as mean (ADC is already summed)
            newcluster.meanposition /= newcluster.sumADC;  // divide by ADC sum to get weighted mean
            newcluster.meanfiber /= newcluster.sumADC;
            newcluster.meantime /= newcluster.NrHits;    // divide by NrHits to get mean
#ifdef HODOSCOPE_FULLTIMESTAMP
            newcluster.meanfulltime /= newcluster.NrHits;    // divide by NrHits to get mean
#endif
            // keep pointers to maximum adc cluster per plane updated
            if (MaxCluster[newcluster.plane-1] == -1 ||
                  fClusterbank[MaxCluster[newcluster.plane-1]].sumADC<newcluster.sumADC)
               MaxCluster[newcluster.plane-1] = fClusterbank.size();

            fClusterbank.push_back(newcluster);
            NrClustersPerPlane[newcluster.plane-1]++;
         }
      } while (!done);
   }

   if (fParam->Debug) PrintAllClusters();

   // Fill eventvise histograms
   if (maxFiber[0]>0 && maxFiber[1]>0) {
      hMaxCh->Fill(maxCh[0]); hMaxCh->Fill(maxCh[1]);
      hMaxFiberall->Fill(maxFiber[0]); hMaxFiberall->Fill(maxFiber[1]+64);
      hMaxFiberpl1->Fill(maxFiber[0]);
      hMaxFiberpl2->Fill(maxFiber[1]);
      hADC->Fill(maxADC[0]); hADC->Fill(maxADC[1]);
      hMaxADC->Fill(maxCh[0],maxADC[0]); hMaxADC->Fill(maxCh[1],maxADC[1]);
      hMaxADCraw->Fill(maxCh[0],maxADCraw[0]); hMaxADCraw->Fill(maxCh[1],maxADCraw[1]);
      hMaxFiber2D->Fill(maxFiber[0],maxFiber[1]);
   }

   hNrClusters->Fill(fClusterbank.size());

   hNrClusters2D->Fill(NrClustersPerPlane[0], NrClustersPerPlane[1]);

   // TODO: for the moment only single hit is detected based on maximum ADC
   // later one could try to extract several positions

   if (MaxCluster[0]>=0 && MaxCluster[1]>=0) {

      hMaxClusterTimediff->Fill(fClusterbank[MaxCluster[0]].meantime-
            fClusterbank[MaxCluster[1]].meantime);

      if (TMath::Abs(fClusterbank[MaxCluster[0]].meantime-
            fClusterbank[MaxCluster[1]].meantime) < fParam->MaxTimeDiff) {

         TFiberHodEvent::hit_t hh;

         hh.X = fClusterbank[MaxCluster[0]].meanposition;
         hh.Y = fClusterbank[MaxCluster[1]].meanposition;

         hh.time = (fClusterbank[MaxCluster[0]].meantime +
                     fClusterbank[MaxCluster[1]].meantime) / 2;

#ifdef HODOSCOPE_FULLTIMESTAMP
         hh.fulltime = (fClusterbank[MaxCluster[0]].meanfulltime +
                            fClusterbank[MaxCluster[1]].meanfulltime) / 2;
#endif
         hh.NumHits = fClusterbank[MaxCluster[0]].NrHits +
                      fClusterbank[MaxCluster[1]].NrHits;

         hh.ADC = fClusterbank[MaxCluster[0]].sumADC +
                  fClusterbank[MaxCluster[1]].sumADC;



         hPosition2D->Fill(hh.X, hh.Y);

         fHodEvent->fHits.push_back(hh);
      }
   }

   // Single Event Display:
   // Fill cumulative single event hit- and cluster histograms only if activated
   if (fParam->SingleEventMode) {
      fSingleEventSlice++;
      if (fSingleEventSlice>100) fSingleEventSlice=1;

      if (hSingleEventCluster2D!=0) {
         // Fill cumulative single hit histogram
         for (Int_t i=0; i<=129; i++)  // include under- and overflow bin
            hSingleEventHit2D->SetBinContent(fSingleEventSlice,i,0);
         for (unsigned ihit=0; ihit<fHitbank.size(); ihit++)
            hSingleEventHit2D->SetBinContent(fSingleEventSlice,
                  fHitbank[ihit].fiber+64*(fHitbank[ihit].plane-1),fHitbank[ihit].ADC);
      }

      // Fill cumulative single cluster histograms
      for (unsigned icluster=1; icluster<=5; icluster++) { // up to 5 clusters are displayed

         if (hSingleEventCluster2D[icluster-1]==0) continue;

         for (Int_t i=0; i<=129; i++)
            hSingleEventCluster2D[icluster-1]->SetBinContent(fSingleEventSlice,i,0);

         if ((icluster-1) < fClusterbank.size()) {
/*            for (unsigned ihit=0; ihit<fClusterbank[icluster-1].hitlist.size(); ihit++) {
               hSingleEventCluster2D[icluster-1]->SetBinContent(fSingleEventSlice,
                     fClusterbank[icluster-1].hitlist[ihit]->fiber+
                     64*(fClusterbank[icluster-1].plane-1),
                     fClusterbank[icluster-1].hitlist[ihit]->ADC);
            }
*/
            hSingleEventCluster2D[icluster-1]->SetBinContent(fSingleEventSlice,0,
                  fClusterbank[icluster-1].meanfiber+
                  64*(fClusterbank[icluster-1].plane-1));
         }
      }
   }
}


void TFiberHodProc::PrintAllHits()
{
   printf("\nDump of Hitbank with %d hits \n", (int) fHitbank.size());
   for (unsigned ihit=0; ihit<fHitbank.size(); ihit++)
      printf("hit %u: fiber %d, plane %d, ADC %f, position:%f, Status %d \n",ihit+1,
            fHitbank[ihit].fiber, fHitbank[ihit].plane, fHitbank[ihit].ADC, fHitbank[ihit].position, fHitbank[ihit].status);
}

void TFiberHodProc::PrintAllClusters()
{
   printf("\nDump of clusterbank with %d clusters\n", (int) fClusterbank.size());
   for (unsigned icluster=0; icluster<fClusterbank.size(); icluster++) {
      printf("Cluster %d: ADCsum:%6.1f, TDCmean %6.1f, Position %6.1f, NrHits: %d \n",
            icluster,fClusterbank[icluster].sumADC, fClusterbank[icluster].meantime, fClusterbank[icluster].meanposition, fClusterbank[icluster].NrHits);
   }
}