TruthHelp.cxx

Go to the documentation of this file.

#include <cassert>
#include <vector>

#include "TMath.h"

#include "AtNuEvent/AtmosEvent.h"
#include "AtNuEvent/AtmosMC.h"
#include "AtNuEvent/AtmosScintHit.h"
#include "AtNuEvent/AtmosStrip.h"
#include "AtNuEvent/AtmosShower.h"
#include "AtNuEvent/AtmosTrack.h"

#include "MessageService/MsgService.h"

#include "AtNuUtils/UtilMisc.h"
#include "AtNuUtils/TruthHelp.h"

class SimpleScintHit
{
 public:
  SimpleScintHit() {TCInd = 0; Strip = 0;}
  SimpleScintHit(int i, int s) {TCInd = i; Strip = s;}
  int TCInd;
  int Strip;
};

class SimpleScintPlane
{
 public:
  SimpleScintPlane() {Plane = 0; TotalDE = 0.0;}
  SimpleScintPlane(int pln) {Plane = pln; TotalDE = 0.0;}
  int Plane;
  double TotalDE;

  vector<SimpleScintHit> Hits;
};

CVSID("$Id: TruthHelp.cxx,v 1.4 2008/04/04 21:52:53 bspeak Exp $");

TruthHelp::TruthHelp()
{
  this->Zero();
}

TruthHelp::TruthHelp(const AtmosEvent* event)
{
  this->Zero();
  this->Fill(event);
}

void TruthHelp::Zero() {
  TrkIdToId.clear();
  AssociateHits.clear();

  AllHits.clear();
  LeptonHits.clear();
  HadronHits.clear();
  NoiseHits.clear();
  DeMuxFailures.clear();
}

void TruthHelp::Fill(const AtmosEvent* event) {
  if(!event) return;
  this->Zero();
  MSG("FarDetEvent",Msg::kVerbose) << "TruthHelp::Fill" << endl;

  MakeTrkIdToId(event);

  MakeAssociateHits(event);

  ClassHits(event);
}

void TruthHelp::MakeTrkIdToId(const AtmosEvent *event) {
  if(!event) return;
  TrkIdToId.clear();
  
  for (int i=0; i<event->ScintHitList->GetEntries(); i++) {
    const AtmosScintHit* hit =
      dynamic_cast<const AtmosScintHit*>(event->ScintHitList->At(i));
    assert(hit);
    //The first hit in the TClonesArray decides if Lepton or not
    if (TrkIdToId.find(TMath::Abs(hit->TrkId)) == TrkIdToId.end()) {
      TrkIdToId[TMath::Abs(hit->TrkId)] = hit->Id;
    }
  }
}

void TruthHelp::MakeAssociateHits(const AtmosEvent *event) {
  if(!event) return;
  AssociateHits.clear();
  DeMuxFailures.clear();

  static SimpleScintPlane ScintPlanes[485];
  for(int i=0; i<485; i++) {
    ScintPlanes[i].Plane = i+1;
    ScintPlanes[i].TotalDE = 0.0;
    ScintPlanes[i].Hits.clear();
  }

  for (int i=0; i<event->ScintHitList->GetEntries(); i++) {
    const AtmosScintHit* hit =
      dynamic_cast<const AtmosScintHit*>(event->ScintHitList->At(i));
    assert(hit);
    ScintPlanes[hit->Plane-1].Hits.push_back(SimpleScintHit(i, hit->Strip));
    ScintPlanes[hit->Plane-1].TotalDE += hit->DE;
  }

  for (int istrip=0; istrip<event->StripList->GetEntries(); istrip++) {
    vector<int> Hit;
    const AtmosStrip *strip = 
      dynamic_cast<AtmosStrip*>(event->StripList->At(istrip));
    SimpleScintPlane &thisScintPlane = ScintPlanes[strip->Plane - 1];

    //Not enough DE in plane noise hit
    if (thisScintPlane.TotalDE < 1e-4) {
      MSG("FarDetEvent",Msg::kDebug) << "Strip " << strip->Strip << " in plane "
             << strip->Plane << " has too little plane DE" << endl;
      AssociateHits.push_back(Hit);
      continue;
    }

    //Check all ganged strips, starting with the bang on strip
    vector<int> Gang0 = UtilStrip::GangedStrips(strip->Strip,0);
    vector<int> Gang1 = UtilStrip::GangedStrips(strip->Strip,1);

    bool MatchActual = false;//True if the actual strip has a match
    for (unsigned int i=0; i<thisScintPlane.Hits.size(); i++) {
      SimpleScintHit &thisScintHit = thisScintPlane.Hits[i];
      vector<int>::iterator F0 = 
        std::find(Gang0.begin(),Gang0.end(),thisScintHit.Strip);
      vector<int>::iterator F1 = 
        std::find(Gang1.begin(),Gang1.end(),thisScintHit.Strip);

      bool MatchStrip = false;
      if(strip->Qadc[0]>0.0 && F0!=Gang0.end() ) MatchStrip = true;
      if(strip->Qadc[1]>0.0 && F1!=Gang1.end() ) MatchStrip = true;
      if(MatchStrip) Hit.push_back(thisScintHit.TCInd);

      if(strip->Qadc[0]>0.0 && F0==Gang0.begin() ) MatchActual = true;
      if(strip->Qadc[1]>0.0 && F1==Gang1.begin() ) MatchActual = true;

      vector<int>::iterator sps;
      if (strip->Qadc[0]>0.0 && F0!=Gang0.begin() && F0!=Gang0.end()) {
        MSG("FarDetEvent",Msg::kSynopsis) << "East DeMux Failure on Plane = " << strip->Plane
               << ", Strip = " << strip->Strip << endl;
        MSG("FarDetEvent",Msg::kSynopsis) << "Ganged Strips =";
        for(unsigned int i=0; i<Gang0.size(); i++) MSG("FarDetEvent",Msg::kSynopsis) << " " << Gang0[i];
        MSG("FarDetEvent",Msg::kSynopsis) << endl;
        MSG("FarDetEvent",Msg::kSynopsis) << "ScintHit Strip = " << thisScintHit.Strip << endl;
      }

      if (strip->Qadc[1]>0.0 && F1!=Gang1.begin() && F1!=Gang1.end()) {
        MSG("FarDetEvent",Msg::kSynopsis) << "West DeMux Failure on Plane = " << strip->Plane
               << ", Strip = " << strip->Strip << endl;
        MSG("FarDetEvent",Msg::kSynopsis) << "Ganged Strips =";
        for(unsigned int i=0; i<Gang1.size(); i++) MSG("FarDetEvent",Msg::kSynopsis) << " " << Gang1[i];
        MSG("FarDetEvent",Msg::kSynopsis) << endl;
        MSG("FarDetEvent",Msg::kSynopsis) << "ScintHit Strip = " << thisScintHit.Strip << endl;
      }
    }
    if(!MatchActual) DeMuxFailures.insert(SetStrip(istrip,strip));
    AssociateHits.push_back(Hit);
  }

  /*
  set<SetStrip>::iterator sitr;
  MSG("FarDetEvent",Msg::kSynopsis) << "DeMuxFailures:" << endl;
  for (sitr=DeMuxFailures.begin(); sitr!=DeMuxFailures.end(); sitr++) {
    MSG("FarDetEvent",Msg::kSynopsis) << "  TCIndex = " << sitr->TCIndex << endl;
  }
  */
}

void TruthHelp::ClassHits(const AtmosEvent *event) {
  if(!event) return;
  AllHits.clear();
  LeptonHits.clear();
  HadronHits.clear();
  NoiseHits.clear();
  
  for (unsigned int i=0; i<AssociateHits.size(); i++) {
    const AtmosStrip* strip =
      dynamic_cast<const AtmosStrip*>(event->StripList->At(i));
    assert(strip);

    AllHits.insert(SetStrip(i, strip));

    const double stripPE = (strip->QPE[0]+strip->QPE[1]);

    //noise hit
    if (AssociateHits[i].size() == 0) {
      NoiseHits.insert(SetStrip(i,stripPE,strip->Ndigits-1));
      continue;
    }

    bool IsLep = false;
    bool IsHad = false;
    
    for (unsigned int j=0; j<AssociateHits[i].size(); j++) {
      int ThisIndex = AssociateHits[i][j];
      const AtmosScintHit* hit =
        dynamic_cast<const AtmosScintHit*>
        (event->ScintHitList->At(ThisIndex));
      assert(hit);
      if(UtilMisc::IdIsLepton(TrkIdToId[TMath::Abs(hit->TrkId)]))
        IsLep = true;
      else IsHad = true;
    }
    if(IsLep) LeptonHits.insert(SetStrip(i,stripPE,strip->Ndigits-1));
    if(IsHad) HadronHits.insert(SetStrip(i,stripPE,strip->Ndigits-1));
  }
}

---