TP0DECALReconModule.cxx

Go to the documentation of this file.

#include "TDatum.hxx"
#include "TGeomInfo.hxx"
#include "TPrincipal.h"
#include "TRealDatum.hxx"
#include "TReconHit.hxx"
#include "TShowerTruthInfo.hxx"
#include "TrackTruthInfo.hxx"

#include "TP0DECALReconModule.hxx"

ClassImp(ND::TP0DECALReconModule::TP0DECALReconTrack);
ClassImp(ND::TP0DECALReconModule::TP0DECALReconShower);

#include "../cvstags/TP0DECALReconModule.cxx"

int GetStackNo(TString Name) {
  if (Name.Contains("LS")) {
    return 1;
  } else if (Name.Contains("LT")) {
    return 2;
  } else if (Name.Contains("LB")) {
    return 3;
  } else if (Name.Contains("RS")) {
    return 4;
  } else if (Name.Contains("RT")) {
    return 5;
  } else if (Name.Contains("RB")) {
    return 6;
  } else {
    return -1;
  }
}

ND::TP0DECALReconModule::TP0DECALReconModule(const char *name,
                                             const char *title) {
  SetNameTitle(name, title);
  // Enable this module by default:
  fIsEnabled = kTRUE;
  fDescription = "P0D ECAL recon output";
  fCVSTagName = CVSTAG;
  fCVSID = CVSID;

  this->SetSplitLevel(1);  // set splitlevel of trees

  // Tree variables
  fReconTracks =
      new TClonesArray("ND::TP0DECALReconModule::TP0DECALReconTrack", 200);
  fReconShowers =
      new TClonesArray("ND::TP0DECALReconModule::TP0DECALReconShower", 200);
  fLayerCharge =
      new TClonesArray("ND::TP0DECALReconModule::TP0DECALLayerCharge", 200);
}

ND::TP0DECALReconModule::~TP0DECALReconModule() {}

Bool_t ND::TP0DECALReconModule::ProcessFirstEvent(ND::TND280Event &event) {
  return true;
}

void ND::TP0DECALReconModule::InitializeModule() {}

void ND::TP0DECALReconModule::InitializeBranches() {
  // Event Level Recon Info
  fOutputTree->Branch("NReconTracks", &fNReconTracks, "NReconTracks/I",fBufferSize)
  ->SetTitle("Number TP0DECALReconTrack Objects ");
  fOutputTree->Branch("NReconShowers", &fNReconShowers, "NReconShowers/I",fBufferSize)
  ->SetTitle("Number of TP0DECALReconShower Objects ");
  fOutputTree->Branch("NLayerCharge", &fNLayerCharge, "NLayerCharge/I",fBufferSize)
  ->SetTitle(" Number of TP0DECALLayerCharge objects");

  // A TClonesArray of TP0DECALReconTrack objects
  fOutputTree->Branch("ReconTracks", &fReconTracks, fBufferSize, fSplitLevel)
  ->SetTitle(" TP0DECALReconTrack Objects");
  // A TClonesArray of TP0DECALReconShower objects
  fOutputTree->Branch("ReconShowers", &fReconShowers, fBufferSize, fSplitLevel)
  ->SetTitle("TP0DECALReconShower Objects ");
  fOutputTree->Branch("LayerCharge", &fLayerCharge, fBufferSize, fSplitLevel)
  ->SetTitle("TP0DECALLayerCharge Objects ");
}

bool ND::TP0DECALReconModule::FillTree(ND::TND280Event &event) {
  bool isMC = event.GetContext().IsMC();

  ///////////////// Add Recon Data ///////////////////////////////////////
  // reset the variables
  fNReconTracks = 0;
  fReconTracks->Clear();
  fNReconShowers = 0;
  fReconShowers->Clear();
  fNLayerCharge = 0;
  fLayerCharge->Clear();

  if (!event.GetFit("P0DECALPid")) {
    // No P0DECal Info - return blank info
    return true;
  }

  ND::THandle<ND::TAlgorithmResult> ecalFinalResult =
      event.GetFit("P0DECALPid");

  for (ND::TDataVector::iterator dataIter = ecalFinalResult->begin();
       dataIter != ecalFinalResult->end(); ++dataIter) {
    std::string class_name = (*dataIter)->ClassName();
    if (class_name.find("ND::TReconObjectContainer") != std::string::npos) {
      ND::TReconObjectContainer *reconContainer =
          dynamic_cast<ND::TReconObjectContainer *>(*dataIter);

      if (!reconContainer) {
        // No P0DECal Info - return blank info
        return true;
      }

      // Should be called final
      TString reconContainerName = reconContainer->GetFullName();

      if (reconContainerName.Contains("final")) {
        // Loop over the top level Recon Container - ReconPid's
        for (ND::TReconObjectContainer::iterator reconIter =
                 reconContainer->begin();
             reconIter != reconContainer->end(); ++reconIter) {
          std::string reco_class_name = (*reconIter)->ClassName();

          // Get the TReconPid from the container
          if (reco_class_name.find("ND::TReconPID") != std::string::npos) {
            ND::THandle<ND::TReconPID> pid = (*reconIter);
            ND::THandle<ND::TReconObjectContainer> conContainer =
                pid->GetConstituents();

            for (ND::TReconObjectContainer::iterator conIter =
                     conContainer->begin();
                 conIter != conContainer->end(); ++conIter) {
              std::string constituent_class_name = (*conIter)->ClassName();

              if (constituent_class_name.find("ND::TReconTrack") !=
                  std::string::npos) {
                ND::THandle<ND::TReconTrack> track = (*conIter);
                TP0DECALReconTrack *ecalReconTrack = 0;

                if (track->GetEDeposit()) {
                  ecalReconTrack =
                      new ((*fReconTracks)[fNReconTracks]) TP0DECALReconTrack;

                  // Fill UniqueID
                  ecalReconTrack->UniqueID = track->GetUniqueID();
                  ecalReconTrack->Direction = track->GetDirection();
                  ecalReconTrack->Position = track->GetPosition();

                  // Module Tagging - TRACKS
                  //-------------------------------------------------------
                  ND::THandle<ND::THitSelection> track_hitsel =
                      track->GetHits();
                  std::string module("unknown");
                  if(track_hitsel){
                    ND::THandle<ND::TReconHit> track_reconhit =
                    *(track_hitsel->begin());
                    module = ND::TGeomInfo::Get()
                    .ECAL()
                    .GetModule(track_reconhit->GetGeomId())
                    .GetName();
                  }
                  ecalReconTrack->Module = module;
                  //-------------------------------------------------------

                  ecalReconTrack->StackNo = (int)track->GetDetectors();

                  // TReconTrack standard Get member functions
                  ecalReconTrack->Curvature = track->GetCurvature();
                  ecalReconTrack->EDeposit = track->GetEDeposit();
                  ecalReconTrack->NDOF = track->GetNDOF();
                  ecalReconTrack->Quality = track->GetQuality();
                  ecalReconTrack->Width = track->GetWidth();

                  // Truth matching info - TRACKS
                  //-----------------------------------------------------
                  int trueID = -999;
                  TLorentzVector trueInitPos(-9999, -9999, -9999, -9999);
                  TLorentzVector trueFinalPos(-9999, -9999, -9999, -9999);
                  TLorentzVector trueInitMom(-9999, -9999, -9999, -9999);
                  int truePDG = 0;
                  int trueParentID = -999;
                  double truePur = -9999.0;
                  double trueEff = -9999.0;

                  ND::THandle<ND::THitSelection> track_hits = track->GetHits();

                  if (isMC) {
                    trueID = TrackTruthInfo::GetG4TrajIDHits(*track_hits,
                                                             trueEff, truePur);
                    ND::THandle<ND::TG4TrajectoryContainer>
                        trajectoryContainer =
                            event.Get<ND::TG4TrajectoryContainer>(
                                "truth/G4Trajectories");
                    ND::THandle<ND::TG4Trajectory> g4traj(0);
                    if (trajectoryContainer) {
                      g4traj = trajectoryContainer->GetTrajectory(trueID);
                    }
                    if (g4traj) {
                      trueInitPos = g4traj->GetInitialPosition();
                      trueFinalPos = g4traj->GetFinalPosition();
                      trueInitMom = g4traj->GetInitialMomentum();
                      truePDG = g4traj->GetPDGEncoding();
                      trueParentID = g4traj->GetParentId();
                    }
                  }

                  ecalReconTrack->TrueID = trueID;
                  ecalReconTrack->TruePDG = truePDG;
                  ecalReconTrack->TrueParentID = trueParentID;
                  ecalReconTrack->TrueHitPurity = truePur;
                  ecalReconTrack->TrueHitEfficiency = trueEff;

                  ecalReconTrack->TrueInitPos = trueInitPos;
                  ecalReconTrack->TrueFinalPos = trueFinalPos;
                  ecalReconTrack->TrueInitMom = trueInitMom;
                  //--------------------------------------------------

                  fNReconTracks++;
                  // Pid Vars
                  if (track->Get<ND::TRealDatum>("AMR")) {
                    ecalReconTrack->AMR =
                        track->Get<ND::TRealDatum>("AMR")->GetValue();
                  } else {
                    ecalReconTrack->AMR = -1;
                  }

                  if (track->Get<ND::TRealDatum>("CWTrackWidth")) {
                    ecalReconTrack->CWTrackWidth =
                        track->Get<ND::TRealDatum>("CWTrackWidth")->GetValue();
                  } else {
                    ecalReconTrack->CWTrackWidth = -1;
                  }

                  if (track->Get<ND::TRealDatum>("Max_Ratio")) {
                    ecalReconTrack->Max_Ratio =
                        track->Get<ND::TRealDatum>("Max_Ratio")->GetValue();
                  } else {
                    ecalReconTrack->Max_Ratio = -1;
                  }

                  if (track->Get<ND::TRealDatum>("NormChargeSD")) {
                    ecalReconTrack->NormChargeSD =
                        track->Get<ND::TRealDatum>("NormChargeSD")->GetValue();
                  } else {
                    ecalReconTrack->NormChargeSD = -1;
                  }

                  if (track->Get<ND::TRealDatum>("PathChargeRatio")) {
                    ecalReconTrack->PathChargeRatio =
                        track->Get<ND::TRealDatum>("PathChargeRatio")
                            ->GetValue();
                  } else {
                    ecalReconTrack->PathChargeRatio = -1;
                  }

                  if (track->Get<ND::TRealDatum>("TrShval")) {
                    ecalReconTrack->TrShval =
                        track->Get<ND::TRealDatum>("TrShval")->GetValue();
                  } else {
                    ecalReconTrack->TrShval = -1;
                  }

                  /// If you need any info from the contributing clusters
                  ND::THandle<ND::TReconObjectContainer> clusterContainer =
                      track->GetConstituents();
                  for (ND::TReconObjectContainer::iterator clusterIter =
                           clusterContainer->begin();
                       clusterIter != clusterContainer->end(); ++clusterIter) {
                    std::string track_constituent_class_name =
                        (*clusterIter)->ClassName();
                    if (track_constituent_class_name.find(
                            "ND::TReconCluster") != std::string::npos) {
                      ND::THandle<ND::TReconCluster> cluster = (*clusterIter);

                      if (cluster->Get<ND::TRealDatum>("BenCWPos")) {

                        ecalReconTrack->BenCWPosition = TVector3(
                          cluster->Get<ND::TRealDatum>("BenCWPos")->at(0),
                          cluster->Get<ND::TRealDatum>("BenCWPos")->at(1),
                          cluster->Get<ND::TRealDatum>("BenCWPos")->at(2)
                          );
                      }

                      // Pid Vars
                      if (cluster->Get<ND::TRealDatum>("AMR")) {
                        ecalReconTrack->AMR =
                            cluster->Get<ND::TRealDatum>("AMR")->GetValue();
                      } else {
                        ecalReconTrack->AMR = -1;
                      }

                      if (cluster->Get<ND::TRealDatum>("CWTrackWidth")) {
                        ecalReconTrack->CWTrackWidth =
                            cluster->Get<ND::TRealDatum>("CWTrackWidth")
                                ->GetValue();
                      } else {
                        ecalReconTrack->CWTrackWidth = -1;
                      }

                      if (cluster->Get<ND::TRealDatum>("Max_Ratio")) {
                        ecalReconTrack->Max_Ratio =
                            cluster->Get<ND::TRealDatum>("Max_Ratio")
                                ->GetValue();
                      } else {
                        ecalReconTrack->Max_Ratio = -1;
                      }

                      if (cluster->Get<ND::TRealDatum>("NormChargeSD")) {
                        ecalReconTrack->NormChargeSD =
                            cluster->Get<ND::TRealDatum>("NormChargeSD")
                                ->GetValue();
                      } else {
                        ecalReconTrack->NormChargeSD = -1;
                      }

                      if (cluster->Get<ND::TRealDatum>("PathChargeRatio")) {
                        ecalReconTrack->PathChargeRatio =
                            cluster->Get<ND::TRealDatum>("PathChargeRatio")
                                ->GetValue();
                      } else {
                        ecalReconTrack->PathChargeRatio = -1;
                      }

                      ecalReconTrack->NHits = (int)cluster->GetHits()->size();
                    }
                  }
                }
              } else if (constituent_class_name.find("ND::TReconShower") !=
                         std::string::npos) {
                ND::THandle<ND::TReconShower> shower = (*conIter);
                TP0DECALReconShower *ecalReconShower = 0;

                if (shower->GetEDeposit()) {
                  ecalReconShower = new ((*fReconShowers)[fNReconShowers])
                      TP0DECALReconShower;

                  ecalReconShower->UniqueID = shower->GetUniqueID();
                  ecalReconShower->StackNo = (int)shower->GetDetectors();

                  // TReconShower standard Get member functions
                  ecalReconShower->ConeAngle = shower->GetConeAngle();
                  ecalReconShower->Direction = shower->GetDirection();
                  ecalReconShower->EDeposit = shower->GetEDeposit();
                  ecalReconShower->NDOF = shower->GetNDOF();
                  ecalReconShower->Position = shower->GetPosition();
                  ecalReconShower->Quality = shower->GetQuality();

                  // Module tagging - SHOWERS
                  //-----------------------------
                  ND::THandle<ND::THitSelection> shower_hitsel =
                      shower->GetHits();
                  std::string module;
                  ND::THandle<ND::TReconHit> shower_reconhit =
                      *(shower_hitsel->begin());
                  module = ND::TGeomInfo::Get()
                               .ECAL()
                               .GetModule(shower_reconhit->GetGeomId())
                               .GetName();
                  ecalReconShower->Module = module;
                  //-----------------------------
                  // Truth matching info - SHOWERS
                  //-----------------------------
                  int trueID = -999;

                  if (isMC) {
                    ND::TShowerTruthInfo showerTruth(pid);
                    std::vector<ND::TShowerTruthInfo::TruthInfo> results =
                        showerTruth.GetPrimaryTrajectories();
                    if (results.size() > 0) {
                      trueID = results[0].ParticleID;
                    }
                  }
                  ecalReconShower->TrueID = trueID;
                  //-----------------------------

                  fNReconShowers++;

                  // Pid Vars - Get these from the consituents of which there
                  // should be just one for p0decal objects
                  if (shower->Get<ND::TRealDatum>("AMR")) {
                    ecalReconShower->AMR =
                        shower->Get<ND::TRealDatum>("AMR")->GetValue();
                  } else {
                    ecalReconShower->AMR = -1;
                  }

                  if (shower->Get<ND::TRealDatum>("CWTrackWidth")) {
                    ecalReconShower->CWTrackWidth =
                        shower->Get<ND::TRealDatum>("CWTrackWidth")->GetValue();
                  } else {
                    ecalReconShower->CWTrackWidth = -1;
                  }

                  if (shower->Get<ND::TRealDatum>("Max_Ratio")) {
                    ecalReconShower->Max_Ratio =
                        shower->Get<ND::TRealDatum>("Max_Ratio")->GetValue();
                  } else {
                    ecalReconShower->Max_Ratio = -1;
                  }

                  if (shower->Get<ND::TRealDatum>("NormChargeSD")) {
                    ecalReconShower->NormChargeSD =
                        shower->Get<ND::TRealDatum>("NormChargeSD")->GetValue();
                  } else {
                    ecalReconShower->NormChargeSD = -1;
                  }

                  if (shower->Get<ND::TRealDatum>("PathChargeRatio")) {
                    ecalReconShower->PathChargeRatio =
                        shower->Get<ND::TRealDatum>("PathChargeRatio")
                            ->GetValue();
                  } else {
                    ecalReconShower->PathChargeRatio = -1;
                  }

                  if (shower->Get<ND::TRealDatum>("TrShval")) {
                    ecalReconShower->TrShval =
                        shower->Get<ND::TRealDatum>("TrShval")->GetValue();
                  } else {
                    ecalReconShower->TrShval = -1;
                  }

                  /// If you need any info from the contributing clusters
                  ND::THandle<ND::TReconObjectContainer> clusterContainer =
                      shower->GetConstituents();
                  for (ND::TReconObjectContainer::iterator clusterIter =
                           clusterContainer->begin();
                       clusterIter != clusterContainer->end(); ++clusterIter) {
                    std::string shower_constituent_class_name =
                        (*clusterIter)->ClassName();
                    if (shower_constituent_class_name.find(
                            "ND::TReconCluster") != std::string::npos) {
                      ND::THandle<ND::TReconCluster> cluster = (*clusterIter);


                      if (cluster->Get<ND::TRealDatum>("BenCWPos")) {
                        ecalReconShower->BenCWPosition = TVector3(
                          cluster->Get<ND::TRealDatum>("BenCWPos")->at(0),
                          cluster->Get<ND::TRealDatum>("BenCWPos")->at(1),
                          cluster->Get<ND::TRealDatum>("BenCWPos")->at(2)
                          );
                      }

                      // Pid Vars -Get these from the consituents of which there
                      // should be just one for p0decal objects
                      if (cluster->Get<ND::TRealDatum>("AMR")) {
                        ecalReconShower->AMR =
                            cluster->Get<ND::TRealDatum>("AMR")->GetValue();
                      } else {
                        ecalReconShower->AMR = -1;
                      }

                      if (cluster->Get<ND::TRealDatum>("CWTrackWidth")) {
                        ecalReconShower->CWTrackWidth =
                            cluster->Get<ND::TRealDatum>("CWTrackWidth")
                                ->GetValue();
                      } else {
                        ecalReconShower->CWTrackWidth = -1;
                      }

                      if (cluster->Get<ND::TRealDatum>("Max_Ratio")) {
                        ecalReconShower->Max_Ratio =
                            cluster->Get<ND::TRealDatum>("Max_Ratio")
                                ->GetValue();
                      } else {
                        ecalReconShower->Max_Ratio = -1;
                      }

                      if (cluster->Get<ND::TRealDatum>("NormChargeSD")) {
                        ecalReconShower->NormChargeSD =
                            cluster->Get<ND::TRealDatum>("NormChargeSD")
                                ->GetValue();
                      } else {
                        ecalReconShower->NormChargeSD = -1;
                      }

                      if (cluster->Get<ND::TRealDatum>("PathChargeRatio")) {
                        ecalReconShower->PathChargeRatio =
                            cluster->Get<ND::TRealDatum>("PathChargeRatio")
                                ->GetValue();
                      } else {
                        ecalReconShower->PathChargeRatio = -1;
                      }

                      ecalReconShower->NHits = (int)cluster->GetHits()->size();
                    }
                  }
                }
              }
            }
          }
        }
      } else if (reconContainerName.Contains("charge_by_layer")) {
        for (ND::TReconObjectContainer::iterator reconIter =
                 reconContainer->begin();
             reconIter != reconContainer->end(); ++reconIter) {
          ND::THandle<ND::TReconCluster> cluster = (*reconIter);
          TP0DECALLayerCharge *LayerCharge = 0;

          LayerCharge =
              new ((*fLayerCharge)[fNLayerCharge]) TP0DECALLayerCharge;
          fNLayerCharge++;

          LayerCharge->module_name = "";
          TString algoname = cluster->GetAlgorithmName();
          if (algoname.Contains("Left")) {
            LayerCharge->module_name += "Left";
          }
          if (algoname.Contains("Right")) {
            LayerCharge->module_name += "Right";
          }
          LayerCharge->module_name += " ";

          if (algoname.Contains("Top")) {
            LayerCharge->module_name += "Top";
          }
          if (algoname.Contains("Bottom")) {
            LayerCharge->module_name += "Bottom";
          }
          if (algoname.Contains("Side")) {
            LayerCharge->module_name += "Side";
          }

          if (cluster->Get<ND::TRealDatum>("first_time")) {
            LayerCharge->FirstTime =
                cluster->Get<ND::TRealDatum>("first_time")->GetValue();
          } else {
            LayerCharge->FirstTime = -1;
          }

          if (cluster->Get<ND::TRealDatum>("last_time")) {
            LayerCharge->LastTime =
                cluster->Get<ND::TRealDatum>("last_time")->GetValue();
          } else {
            LayerCharge->LastTime = -1;
          }

          if (cluster->Get<ND::TRealDatum>("layer_charge")) {
            for (int nlay = 0; nlay < 6; nlay++) {
              LayerCharge->Charge[nlay] =
                  cluster->Get<ND::TRealDatum>("layer_charge")->at(nlay);
            }
          }
        }
      }
    }
  }
  return true;
}

ND::TP0DECALReconModule::TP0DECALReconTrack::~TP0DECALReconTrack() {}
ND::TP0DECALReconModule::TP0DECALReconShower::~TP0DECALReconShower() {}
ND::TP0DECALReconModule::TP0DECALLayerCharge::~TP0DECALLayerCharge() {}