MMultipletValidation.cc

Go to the documentation of this file.

// this will be a module that checks for coherence in coincidence multiplets
// and selects just the ones that exist, symmetric, after a given set of cuts
// meaning that both MAIN <-> COINCIDENT pass the cuts
 
#include "MMultipletValidation.hh"
#include "QEvent.hh"
#include "QRawEvent.hh"
#include "QBaseType.hh"
#include "QCoincidenceData.hh"
#include "QBool.hh"
 
#include <string>
#include <map>
#include <set>
 
using namespace std;
using namespace Diana;
 
REGISTER_MODULE(MMultipletValidation);
 
LightMultiplet::LightMultiplet( ULong64_t _tag ){
  fValid = true; // if no-one comes to set this multiplet invalid, assume it OK
  fTag = _tag;
}
 
bool LightMultiplet::Insert( ULong64_t _tag ){
  if( fCoincidentTag.count( _tag ) != 0 )
    return false;
  fCoincidentTag.insert( _tag );
  return true;
}
 
void LightMultiplet::SetInvalid(){
  fValid = false;
}
 
bool LightMultiplet::GetIsValid(){
  return fValid;
}
void MMultipletValidation::Init( Diana::QEvent& ev ){
  fCoincidenceLabel = GetString("CoincidenceLabel","FastCoincidence@CoincidenceData");
  // need to read EventNumber and Channel
  ev.Require<QHeader>("DAQ", "Header");
  ev.Require<QPulseInfo>("DAQ", "PulseInfo");
  ev.RequireByLabel<QCoincidenceData>( fCoincidenceLabel );
  ev.Add<QBool>("Valid");
  ev.Add<QBool>("Selected");
}
 
void MMultipletValidation::Do( Diana::QEvent& ev ){
  const QHeader& header            = ev.Get<QHeader>("DAQ", "Header");
  const QPulseInfo& pulseinfo      = ev.Get<QPulseInfo>("DAQ", "PulseInfo");
  const QCoincidenceData& coindata = ev.GetByLabel<QCoincidenceData>( fCoincidenceLabel );
  const int channel                = pulseinfo.GetChannelId();
  const int tower                  = ( (channel-1)/52 ) + 1;
  const int EventNumber            = header.GetEventNumber();
  ULong64_t tag                    = 100*EventNumber + tower; // unique tag for event
  // 1st iteration: build "graph" in RAM of multiplets
  if( GetIteration() == 1 ){
    if( fMultiplets.count( tag ) != 0 )
      Panic("Something impossible happened: unique tag %llu found twice!",tag);
    fMultiplets[ tag ] = new LightMultiplet( tag );
    // loop over coincident event data
    for(unsigned int j=0; j < coindata.fCoincidentChannels.size(); j++){
      QCoincidentChannel IT = coindata.fCoincidentChannels[j];
      int cointower         = ((IT.fChannelId - 1)/52) + 1;
      ULong64_t thistag     = 100*IT.fEventNumber + cointower;
      if( !(fMultiplets.at(tag))->Insert( thistag ) )
        Panic("Current tag: %llu. Found duplicate tag %llu in coincident tag-list.",tag,thistag);
    }
  }else{// 2nd iteration: perform selection
    ev.Get<QBool>("Valid")  = (fMultiplets.at( tag ))->GetIsValid();
    ev.Get<QBool>("Selected") = true;
  }
}
 
void MMultipletValidation::Done( ){
  if( GetIteration() == 1 ){
    SetRunAgain(true);
    // 1st part: assuming the coincidences are built correctly
    //           and mismatches are due just to different selections
    /* perform all kind of checks and prepare bool flags for next iteration
     * 1. loop over the primary keys (both valid and invalid).
     *    if valid..
     *    a. check that each coincident key exists as primary (was selected)
     *    b. if yes, go ahead with next coincident key.
     *    c. if no, mark the main key as invalid and go to 2.
     * 2. if invalid, loop over its coincidents and set them all invalid
     */
    std::map< ULong64_t, LightMultiplet* >::iterator mainIT;
    for( mainIT = fMultiplets.begin(); mainIT != fMultiplets.end(); ++mainIT )
      {// loop over multiplets (main event tags)
        if( (mainIT->second)->GetIsValid() ){
          std::set<ULong64_t>::iterator coinIT;
          for(coinIT = ((mainIT->second)->fCoincidentTag).begin(); coinIT != ((mainIT->second)->fCoincidentTag).end(); ++coinIT )
            {// loop over coincident keys
              if( fMultiplets.count( *coinIT ) <= 0 ){// coincident key not found as main
                (mainIT->second)->SetInvalid();// set main invalid
                break;// stop checking other coincidents
              }
            }
        }
        // if main is invalid at this point, set invalid all the existing coincidents
        if( false == (mainIT->second)->GetIsValid() ){
          std::set<ULong64_t>::iterator coinIT;
          std::map<ULong64_t, LightMultiplet*>::iterator mainIT2;
          for(coinIT = ((mainIT->second)->fCoincidentTag).begin(); coinIT != ((mainIT->second)->fCoincidentTag).end(); ++coinIT )
            {// loop over coincident keys
              mainIT2 = fMultiplets.find( *coinIT );
              if( mainIT2 != fMultiplets.end() ){// if exists in the main key list...
                (mainIT2->second)->SetInvalid();// ... set it invalid
              }
            }
        }
      }
    // 2nd part: check coherence and remove iteratively mismatched multiplets
    bool myRunAgain;
    int  myRunCounter = 0;
    do{
      myRunCounter++;
      Info("Running 2nd level check of coincident multiplets. Iteration %d", myRunCounter);
      if( myRunCounter > 1 )
        Warn("The 2nd level check had to modify multiplet validity at least once. This means multiplets were NOT symmetric in the first place.");
      myRunAgain = false;
      for( mainIT = fMultiplets.begin(); mainIT != fMultiplets.end(); ++mainIT )
        {// loop over multiplets (main event tags)
          if( (mainIT->second)->GetIsValid() )
            { // if main valid go ahead
              std::set<ULong64_t>::iterator coinIT;
              std::map<ULong64_t, LightMultiplet*>::iterator mainIT2;
              for(coinIT = ((mainIT->second)->fCoincidentTag).begin(); coinIT != ((mainIT->second)->fCoincidentTag).end(); ++coinIT )
                {// loop over coincident keys...
                  if( fMultiplets.count( *coinIT ) <= 0 )
                    Panic("Event tag %llu, coincident of valid main event tag %llu not found.", *coinIT, mainIT->first);
                  mainIT2 = fMultiplets.find( *coinIT ); // find it in the main event list...
                  if( (mainIT2->second)->GetIsValid() && ((mainIT2->second)->fCoincidentTag).count( mainIT->first ) )
                    continue; // if from there you can come back to the main, everything is fine.
                  (mainIT->second)->SetInvalid(); // otherwise the main is invalid ... 
                  if( !myRunAgain )
                    myRunAgain = true; // run another iteration later
                  Info("Event tag %llu found as coincident of valid main event tag %llu, either is invalid or does not list the main tag among its coincidents. Set main invalid.", mainIT2->first, mainIT->first);
                  Info("Setting invalid all coincident events of main event tag %llu", mainIT->first);
                  --mainIT; // and also its coincidents! (this loops back on the same event -> finds it invalid -> invalidates all the coincidents
                  break;
                }
            }else{ // if main invalid check that ALL coincidents are invalid!
            std::set<ULong64_t>::iterator coinIT;
            std::map<ULong64_t, LightMultiplet*>::iterator mainIT2;
            for(coinIT = ((mainIT->second)->fCoincidentTag).begin(); coinIT != ((mainIT->second)->fCoincidentTag).end(); ++coinIT )
              {// loop over coincidents
                mainIT2 = fMultiplets.find( *coinIT );
                if( mainIT2 != fMultiplets.end() ){ // if the coincident exists in the main key list...
                  if( (mainIT2->second)->GetIsValid() ){ // ... and is valid ...
                    (mainIT2->second)->SetInvalid(); // ... set it invalid!
                    if( !myRunAgain )
                      myRunAgain = true; // run another iteration later
                    Info("Event tag %llu found as coincident of invalid main event tag %llu, was set valid. Set invalid.", mainIT2->first, mainIT->first); // and inform the user.
                  }
                }
              }
          }
        }
    }while( myRunAgain );
  }else{
    SetRunAgain(false);
    std::map< ULong64_t, LightMultiplet*>::iterator it;
    for( it = fMultiplets.begin(); it != fMultiplets.end(); ++it)
      delete it->second;
  }
}