KVEventFiltering.cpp

//Created by KVClassFactory on Thu Jun 21 17:01:26 2012
//Author: John Frankland,,,
 
#include "KVEventFiltering.h"
#include "KVMultiDetArray.h"
#include "KV2Body.h"
#include "KVDBSystem.h"
#include "KVDBRun.h"
#include "KVDataSet.h"
#include "KVDataSetManager.h"
#include "KVGeoNavigator.h"
#include "TRandom.h"
 
ClassImp(KVEventFiltering)
 
 
 
 
 
KVEventFiltering::KVEventFiltering()
{
   // Default constructor
   fTransformKinematics = kTRUE;
   fNewFrame = "";
   fRotate = kTRUE;
#ifdef WITH_GEMINI
   fGemini = kFALSE;
   fGemDecayPerEvent = 1;
   fGemAddRotEner = kFALSE;
#endif
#ifdef DEBUG_FILTER
   memory_check = KVClassMonitor::GetInstance();
   SetEventsReadInterval(100);
#endif
}
 
 
 
 
 
KVEventFiltering::KVEventFiltering(const KVEventFiltering& obj)  : KVEventSelector()
{
   // Copy constructor
   // This ctor is used to make a copy of an existing object (for example
   // when a method returns an object), and it is always a good idea to
   // implement it.
   // If your class allocates memory in its constructor(s) then it is ESSENTIAL :-)
 
   fTransformKinematics = kTRUE;
   fNewFrame = "";
   obj.Copy(*this);
}
 
 
 
 
KVEventFiltering::~KVEventFiltering()
{
   // Destructor
#ifdef DEBUG_FILTER
   delete memory_check;
#endif
}
 
 
 
 
 
void KVEventFiltering::Copy(TObject& obj) const
{
   // This method copies the current state of 'this' object into 'obj'
   // You should add here any member variables, for example:
   //    (supposing a member variable KVEventFiltering::fToto)
   //    CastedObj.fToto = fToto;
   // or
   //    CastedObj.SetToto( GetToto() );
 
   KVEventSelector::Copy(obj);
   KVEventFiltering& CastedObj = (KVEventFiltering&)obj;
   CastedObj.fRotate = fRotate;
#ifdef WITH_GEMINI
   CastedObj.fGemini = fGemini;
   CastedObj.fGemDecayPerEvent = fGemDecayPerEvent;
#endif
}
 
 
 
 
void KVEventFiltering::RandomRotation(KVEvent* to_rotate, const TString& frame_name) const
{
   // do random phi rotation around z-axis
   // if frame_name is given, apply rotation to that frame
   //
   // store phi rotation angle [radians] in event parameter "RANDOM_PHI"
   TRotation r;
   Double_t phi = gRandom->Uniform(TMath::TwoPi());
   r.RotateZ(phi);
   if (frame_name != "") to_rotate->SetFrame("rotated_frame", frame_name, r);
   else to_rotate->SetFrame("rotated_frame", r);
   to_rotate->SetParameter("RANDOM_PHI", phi);
}
 
 
 
 
Bool_t KVEventFiltering::Analysis()
{
   // Event-by-event filtering of simulated data.
   // If needed (fTransformKinematics = kTRUE), kinematics of event are transformed
   // to laboratory frame using C.M. velocity calculated in InitAnalysis().
   // Detection of particles in event is simulated with KVMultiDetArray::DetectEvent,
   // then the reconstructed detected event is treated by the same identification and calibration
   // procedures as for experimental data.
 
   // few event counter print at the beginning to be sure the process started properly
   // because in case GEMINI decay is used, it sometimes stops (randomly) after few events
   if ((fEVN <= 10)) Info("Analysis", "%d events processed", (int)fEVN);
   else if ((fEVN <= 1000) && !(fEVN % 100)) Info("Analysis", "%d events processed", (int)fEVN);
 
#ifdef DEBUG_FILTER
   if (fEVN == 2500) memory_check->SetInitStatistics();
   else if (fEVN > 4995) memory_check->CompareToInit();
#endif
   KVEvent* to_be_detected = GetEvent();
#ifdef WITH_GEMINI
   Int_t iterations = 1;
   if (fGemini) iterations = fGemDecayPerEvent;
#endif
   if (to_be_detected->GetNumber()) to_be_detected->SetParameter("SIMEVENT_NUMBER", (int)to_be_detected->GetNumber());
   to_be_detected->SetParameter("SIMEVENT_TREE_ENTRY", (int)fTreeEntry);
#ifdef WITH_GEMINI
   do {
      if (fGemini) {
         try {
            GEM.DecayEvent((KVSimEvent*)GetEvent(), &fGemEvent, fGemAddRotEner);
         }
         catch (...) {
            continue;
         }
         //Copy any parameters associated with simulated event into the Gemini-decayed event
         GetEvent()->GetParameters()->Copy(*(fGemEvent.GetParameters()));
         to_be_detected = &fGemEvent;
      }
#endif
      if (fTransformKinematics) {
         if (fNewFrame == "proj")   to_be_detected->SetFrame("lab", fProjVelocity);
         else                    to_be_detected->SetFrame("lab", fCMVelocity);
      }
      else {
         // default kinematics is lab frame. make sure name is set.
         to_be_detected->SetFrameName("lab");
      }
      if (fRotate) {
         RandomRotation(to_be_detected, "lab");
         fDetSimulator.DetectEvent(to_be_detected, "rotated_frame");
      }
      else {
         fDetSimulator.DetectEvent(to_be_detected, "lab");
      }
 
      fEventReconstructor->SetDetectionFrame(fDetSimulator.GetDetectionFrame());
      fEventReconstructor->SetSimEvent(to_be_detected);
      fEventReconstructor->ReconstructEvent();
 
      fReconEvent->SetNumber(fEVN++);
      fReconEvent->SetFrameName("lab");
      fReconEvent->GetParameters()->Concatenate(*to_be_detected->GetParameters());
 
      FillTree();
#ifdef WITH_GEMINI
   }
   while (fGemini && --iterations);
#endif
 
 
   return kTRUE;
}
 
 
 
 
void KVEventFiltering::EndAnalysis()
{
   gEnv->SetValue(Form("%s.HasCalibIdentInfos", GetOpt("Dataset").Data()), fIdCalMode.Data());
}
 
 
 
 
void KVEventFiltering::EndRun()
{
}
 
 
 
 
void KVEventFiltering::InitAnalysis()
{
   // Select required dataset for filtering (option "Dataset")
   // Build the associated multidetector geometry.
   // Calculate C.M. velocity associated with required experimental collision
   // kinematics (option "System").
   //
   // Set the parameters of the detectors according to the required run
   // if given (option "Run"), or the first run of the given system otherwise.
   //
   // Open file for filtered data (see KVEventFiltering::OpenOutputFile), which will
   // be stored in a TTree with name 'ReconstructedEvents', in a branch with name
   // 'ReconEvent'. The class used for reconstructed events depends on the dataset,
   // it is given by KVDataSet::GetReconstructedEventClassName().
 
   if (fDisableCreateTreeFile) return; //when running with PROOF, only execute on workers
 
   TString dataset = GetOpt("Dataset").Data();
   if (!gDataSetManager) {
      gDataSetManager = new KVDataSetManager;
      gDataSetManager->Init();
   }
   gDataSetManager->GetDataSet(dataset)->cd();
 
   TString system;
   if (IsOptGiven("System")) system = GetOpt("System").Data();
   KVDBSystem* sys = (gExpDB && system != "" ? gExpDB->GetSystem(system) : nullptr);
 
   KV2Body* tb = nullptr;
 
   Bool_t justcreated = kFALSE;
   if (sys) tb =  sys->GetKinematics();
   else if (system != "") {
      tb = new KV2Body(system);
      tb->CalculateKinematics();
      justcreated = kTRUE;
   }
 
   fCMVelocity = (tb ? tb->GetCMVelocity() : TVector3(0, 0, 0));
   fCMVelocity *= -1.0;
 
   fProjVelocity = (tb ? tb->GetNucleus(1)->GetVelocity() : TVector3(0, 0, 0));
   fProjVelocity *= -1.0;
 
   if (justcreated)
      delete tb;
 
   Int_t run = 0;
   if (IsOptGiven("Run")) {
      run = GetOpt("Run").Atoi();
      Info("InitAnalysis", "Run given in options = %d", run);
   }
   if (!run) {
      if (sys) {
         run = ((KVDBRun*)sys->GetRuns().First())->GetNumber();
         Info("InitAnalysis", "Using first run for system = %d", run);
      }
      else {
         Info("InitAnalysis", "No run information");
         run = -1;
      }
   }
 
   fIdCalMode = gEnv->GetValue(Form("%s.HasCalibIdentInfos", dataset.Data()), "no");
 
   TString filt = GetOpt("Filter").Data();
   if (filt == "GeoThresh") {
      // in geo+thresholds mode, we ignore any calibration/identification parameters for the dataset,
      // i.e. all detectors and all identification ntelescopes are supposed to work
      gEnv->SetValue(Form("%s.HasCalibIdentInfos", dataset.Data()), "no");
   }
 
   KVMultiDetArray::MakeMultiDetector(dataset, run);
   if (run == -1) gMultiDetArray->InitializeIDTelescopes();
 
   fDetSimulator.SetArray(gMultiDetArray);
 
   if (filt == "Geo") {
      gMultiDetArray->SetFilterType(KVMultiDetArray::kFilterType_Geo);
      Info("InitAnalysis", "Geometric filter");
   }
   else if (filt == "GeoThresh") {
      gMultiDetArray->SetFilterType(KVMultiDetArray::kFilterType_GeoThresh);
      Info("InitAnalysis", "Geometry + thresholds filter");
   }
   else if (filt == "Full") {
      gMultiDetArray->SetFilterType(KVMultiDetArray::kFilterType_GeoThresh);
      Info("InitAnalysis", "Geometry + thresholds filter + implemented identifications and calibrations.");
   }
   TString kine = GetOpt("Kinematics").Data();
   if (kine == "lab") {
      fTransformKinematics = kFALSE;
      Info("InitAnalysis", "Simulation is in laboratory/detector frame");
   }
   else {
      fNewFrame = kine;
      Info("InitAnalysis", "Simulation will be transformed to laboratory/detector frame");
   }
 
   if (IsOptGiven("PhiRot")) {
      if (GetOpt("PhiRot") == "no") fRotate = kFALSE;
   }
   if (fRotate) Info("InitAnalysis", "Random phi rotation around beam axis performed for each event");
#ifdef WITH_GEMINI
   if (IsOptGiven("Gemini")) {
      if (GetOpt("Gemini") == "yes") fGemini = kTRUE;
      fGemAddRotEner = kFALSE;
      if (IsOptGiven("GemDecayPerEvent")) fGemDecayPerEvent = GetOpt("GemDecayPerEvent").Atoi();
      else fGemDecayPerEvent = 1;
      if (IsOptGiven("GemAddRotEner")) {
         if (GetOpt("GemAddRotEner") == "yes") fGemAddRotEner = kTRUE;
      }
   }
   if (fGemini) {
      Info("InitAnalysis", "Statistical decay with Gemini++ for each event before detection");
      if (fGemDecayPerEvent > 1) Info("InitAnalysis", "  -- %d decays per primary event", fGemDecayPerEvent);
      if (fGemAddRotEner) Info("InitAnalysis", "  -- Rotational energy will be added to excitation energy");
   }
#endif
   if (filt == "Full" && !gDataSet->HasCalibIdentInfos()) {
      // for old data without implementation of calibration/identification
      // we set status of identifications according to experimental informations in IdentificationBilan.dat
      // this "turns off" any telescopes which were not working during the experimental run
 
      // beforehand we have to "turn on" all identification telescopes by initializing them
      // as this will not have been done yet
      gMultiDetArray->InitializeIDTelescopes();
      // this will also set informations on (simulated) mass identification capabilities for
      // certain telescopes
 
      TString fullpath;
      if (sys && KVBase::SearchKVFile("IdentificationBilan.dat", fullpath, gDataSet->GetName())) {
         Info("InitAnalysis", "Setting identification bilan...");
         TString sysname = sys->GetBatchName();
         KVIDTelescope::OpenIdentificationBilan(fullpath);
         TIter next(gMultiDetArray->GetListOfIDTelescopes());
         KVIDTelescope* idt;
         while ((idt = (KVIDTelescope*)next())) {
            idt->CheckIdentificationBilan(sysname);
         }
      }
   }
   gMultiDetArray->PrintStatusOfIDTelescopes();
 
   OpenOutputFile(sys, run);
   TTree* t{nullptr};
   if (sys) t = AddTree("ReconstructedEvents", Form("%s filtered with %s (%s)", GetOpt("SimTitle").Data(), gMultiDetArray->GetTitle(), sys->GetName()));
   else t = AddTree("ReconstructedEvents", Form("%s filtered with %s", GetOpt("SimTitle").Data(), gMultiDetArray->GetTitle()));
 
   TString reconevclass = gDataSet->GetReconstructedEventClassName();
   fReconEvent = (KVReconstructedEvent*)TClass::GetClass(reconevclass)->New();
   KVEvent::MakeEventBranch(t, "ReconEvent", fReconEvent);
 
   auto fer = new KVFilterEventReconstructor(gMultiDetArray, fReconEvent);
 
   // look for required data quality audit given as
   //   [dataset][/[system]]
   // [dataset] is not necessarily the same as the one used for detector geometry
   if (IsOptGiven("DataQualityAudit"))
   {
      KVString audit = GetOpt("DataQualityAudit");
      audit.Begin("/");
      KVString dataset,system;
      dataset = audit.Next();
      if(!audit.End())
         system = audit.Next();
      auto dqaf = TFile::Open(KVDataSet::GetFullPathToDataSetFile(dataset,"DataQualityAudit.root"));
      auto dqa = system.IsNull() ? dqaf->Get<KVDataQualityAudit>(dataset) : dqaf->Get<KVDataQualityAudit>(system);
      fer->SetDataQualityAudit(dqa);
 
      if(IsOptGiven("ExtensibleID"))
      {
         // the ID telescopes in the ExtensibleID list can go beyond the largest identified Z in the audit
         // using theoretical (calculated) identification thresholds
         KVString idlabs = GetOpt("ExtensibleID");
         idlabs.Begin("|");
         while(!idlabs.End())
         {
            KVFilterGroupReconstructor::EnableIDExtensionBeyondAudit(idlabs.Next());
         }
      }
   }
 
   fEventReconstructor.reset(fer);
}
 
 
 
 
void KVEventFiltering::InitRun()
{
   fEVN = 0;
#ifdef DEBUG_FILTER
   SetEventsReadInterval(100);
#endif
}
 
 
 
 
void KVEventFiltering::OpenOutputFile(KVDBSystem* S, Int_t run)
{
   // Open ROOT file for new filtered events TTree.
   //
   // The filename is built up from the original simulation filename and the values
   // of various options:
   //
   //       [simfile]_[Gemini]_filt=[filter-type]_[dataset]_[system]_run=[run-number].root
   //
   // In addition, informations on the filtered data are stored in the file as
   // TNamed objects. These can be read by KVSimDir::AnalyseFile:
   //
   // KEY: TNamed System;1 title=[full system name]
   // KEY: TNamed Dataset;1   title=[dataset name]
   // KEY: TNamed Run;1 title=[run-number]
   // KEY: TNamed Filter;1 title=[filter-type]
   // KEY: TNamed Origin;1 title=[name of simulation file]
   // KEY: TNamed RandomPhi;1 title=[yes/no, random rotation about beam axis]
   // KEY: TNamed Gemini++;1 title=[yes/no, Gemini++ decay before detection]
   // KEY: TNamed GemDecayPerEvent;1 title=[number of Gemini++ decays per primary event]
   // KEY: TNamed GemAddRotEner;1 title=[Enable or not the addition of the rotational energy to the excitation energy]
   //
   TString basefile = GetOpt("SimFileName");
   basefile.Remove(basefile.Index(".root"), 5);
   TString outfile = basefile;
#ifdef WITH_GEMINI
   if (fGemini) {
      outfile += "_Gemini";
      if (fGemDecayPerEvent > 1) outfile += fGemDecayPerEvent;
      if (fGemAddRotEner) outfile += "AddedRotEner";
   }
#endif
   outfile += "_filt=";
   outfile += GetOpt("Filter");
   outfile += "_";
   outfile += gDataSet->GetName();
 
   if (S) {
      outfile += "_";
      outfile += S->GetBatchName();
   }
   else if (IsOptGiven("System")) {
      TString tmp = GetOpt("System");
      tmp.ReplaceAll(" ", "");
      tmp.ReplaceAll("@", "");
      tmp.ReplaceAll("MeV/A", "");
      tmp.ReplaceAll("+", "");
      outfile += "_";
      outfile += tmp.Data();
   }
   if (run > 0) {
      outfile += "_run=";
      outfile += Form("%d", run);
   }
   if (IsOptGiven("DataQualityAudit"))
   {
      KVString tmp(GetOpt("DataQualityAudit"));
      tmp.ReplaceAll("/",":");
      outfile += "_DQA=" + tmp;
   }
   outfile += ".root";
 
   TString fullpath;
   AssignAndDelete(fullpath, gSystem->ConcatFileName(GetOpt("OutputDir").Data(), outfile.Data()));
 
   SetJobOutputFileName(fullpath);
 
   TDirectory* curdir = gDirectory;
   writeFile->cd();
   if (S) {
      TNamed* system = new TNamed("System", S->GetName());
      system->Write();
   }
   else if (GetOpt("System"))(new TNamed("System", GetOpt("System").Data()))->Write();
 
   (new TNamed("Dataset", gDataSet->GetName()))->Write();
   if (run > 0)(new TNamed("Run", Form("%d", run)))->Write();
   (new TNamed("Filter", GetOpt("Filter").Data()))->Write();
   if (IsOptGiven("DataQualityAudit"))(new TNamed("DataQualityAudit", GetOpt("DataQualityAudit").Data()))->Write();
   (new TNamed("Origin", (basefile + ".root").Data()))->Write();
   (new TNamed("RandomPhi", (fRotate ? "yes" : "no")))->Write();
#ifdef WITH_GEMINI
   (new TNamed("Gemini++", (fGemini ? "yes" : "no")))->Write();
   (new TNamed("GemDecayPerEvent", Form("%d", fGemDecayPerEvent)))->Write();
   (new TNamed("GemAddRotEner", (fGemAddRotEner ? "yes" : "no")))->Write();
#endif
   curdir->cd();
}