KVINDRADSTConverter.cpp

#include "KVDataAnalyser.h"
#define KVDSTTreeConverter_cxx
 
#include "KVINDRADSTConverter.h"
#include "KVUserAnalysisOptionList.h"
#include "KVINDRADB.h"
#include <KVDetectorStack.h>
#include <KVReconNucTrajectory.h>
 
 
 
void KVINDRADSTConverter::Begin(TTree* /*tree*/)
{
   // The Begin() function is called at the start of the query.
   // When running with PROOF Begin() is only called on the client.
   // The tree argument is deprecated (on PROOF 0 is passed).
 
   Info("Begin", "Reading %lld events from file", total_events);
   events_read = 0;
   KVUserAnalysisOptionList option_list;
   option_list.ParseOptions(GetOption());
   run_number = option_list.GetOpt("run_number").Atoi();
   run_index.set_run(run_number);
   if (!gIndra)
      KVMultiDetArray::MakeMultiDetector(gDataSet->GetName(), run_number);
   else
      gIndra->SetParameters(run_number);
 
   with_phoswich = gIndra->GetListOfPhoswich() && !gIndra->GetListOfPhoswich()->IsEmpty();
   camp1 = gDataSet->IsCalled("INDRA_camp1");
   camp2 = gDataSet->IsCalled("INDRA_camp2");
   camp4 = gDataSet->IsCalled("INDRA_camp4");
 
   if (!the_event) the_event = new KVReconstructedEvent;
 
   output_file = gDataSet->NewRunfile("recon", run_index);
   output_tree = new TTree("ReconEvents", Form("%s : %s : root events converted from DST",
                           gIndraDB->GetDBRunFile(run_index).GetName(),
                           gIndraDB->GetDBRunFile(run_index).GetTitle()));
   KVEvent::MakeEventBranch(output_tree, "ReconEvent", the_event);
}
 
 
 
 
void KVINDRADSTConverter::SlaveBegin(TTree* /*tree*/)
{
   // The SlaveBegin() function is called after the Begin() function.
   // When running with PROOF SlaveBegin() is called on each slave server.
   // The tree argument is deprecated (on PROOF 0 is passed).
 
   TString option = GetOption();
}
 
 
 
 
Bool_t KVINDRADSTConverter::Process(Long64_t entry)
{
   fReader.SetLocalEntry(entry);
 
   the_event->Clear();
 
   for (int i = 0; i < *npart_traite; ++i) {
      if (idcode[i] >= KVINDRA::NO_IDENTIFICATION) continue;
 
      if (idcode[i] == KVINDRA::ID_GAMMA) {
         KVDetector* CSI = nullptr;
         if (with_phoswich && icou[i] == 1)
            CSI = get_phos(i);
         else
            CSI = get_csi(i);
         if (!CSI) {
            Error("Process", "Got no CSI for cod=%d cou=%d mod=%d", idcode[i], icou[i], imod[i]);
         }
         else {
            the_event->GetParameters()->IncrementValue("INDRA_GAMMA_MULT", 1);
            the_event->GetParameters()->IncrementValue("INDRA_GAMMA_DETS", CSI->GetName());
         }
         continue;
      }
      KVINDRADetector* det = nullptr;
 
      if (stopped_in_csi(i)) det = get_csi(i);
      else if (stopped_in_si(i)) det = get_si(i);
      else if (stopped_in_ci(i)) det = get_ci(i);
      else if (stopped_in_phos(i)) det = get_phos(i);
      else if (stopped_in_si75(i)) det = get_si75(i);
      else if (stopped_in_sili(i)) det = get_sili(i);
      if (!det) {
         // nothing fired? try direct approach?
         switch (idcode[i]) {
            case KVINDRA::ID_CSI_PSA:
            case KVINDRA::ID_CSI_FRAGMENT:
               if (!with_phoswich || icou[i] > 1)
                  det = gIndra->GetDetectorByType((UInt_t)icou[i], (UInt_t)imod[i], CsI_T);
               else
                  det = gIndra->GetDetectorByType((UInt_t)icou[i], (UInt_t)imod[i], Phos_T);
               break;
 
            default:
               det = nullptr;
         }
      }
 
      auto idtelescope_code = idcode[i];
      switch (idcode[i]) {
         case KVINDRA::ID_STOPPED_IN_FIRST_STAGE:
            if (!with_phoswich && icou[i] == 1)
               idtelescope_code = KVINDRA::ID_SI_CSI;
            else if (icou[i] < 10)
               idtelescope_code = KVINDRA::ID_CI_SI;
            else if (icou[i] >= 10) {
               if (si75_fired(i)) idtelescope_code = KVINDRA::ID_CI_SI75;
               else idtelescope_code = KVINDRA::ID_CI_CSI;
            }
            break;
 
         case KVINDRA::ID_CI_SI_COHERENCY:
            if (icou[i] > 1 && icou[i] < 10) {
               idtelescope_code = KVINDRA::ID_CI_SI;
            }
            else if (icou[i] > 9) {
               idtelescope_code = KVINDRA::ID_CI_CSI;
               det = get_csi(i);
            }
            break;
 
         case KVINDRA::ID_NEUTRON:
            if (icou[i] < 10) {
               idtelescope_code = KVINDRA::ID_SI_CSI;
            }
            break;
 
         case KVINDRA::ID_CI_MULTIHIT:
            if (icou[i] > 1 && icou[i] < 10) {
               idtelescope_code = KVINDRA::ID_CI_SI;
            }
            else if (icou[i] > 9) {
               idtelescope_code = KVINDRA::ID_CI_CSI;
            }
            break;
 
         case KVINDRA::ID_CI_COHERENCY:
            if (icou[i] > 1 && icou[i] < 10) {
               idtelescope_code = KVINDRA::ID_CI_SI;
               det = get_si(i);
            }
            else if (icou[i] > 9) {
               idtelescope_code = KVINDRA::ID_CI_CSI;
               det = get_csi(i);
            }
            break;
 
         case KVINDRA::ID_CSI_FRAGMENT:
            if (with_phoswich && icou[i] == 1) {
               idtelescope_code = KVINDRA::ID_PHOSWICH;
            }
            else {
               idtelescope_code = KVINDRA::ID_CSI_PSA;
            }
            break;
 
         case KVINDRA::ID_CSI_MASS_OUT_OF_RANGE:
            if (icou[i] > 1) {
               idtelescope_code = KVINDRA::ID_CSI_PSA;
            }
            break;
      }
      if (!det) {
         Warning("Process",
                 "Got no detector with id=%d cou=%d mod=%d z=%d a=%d ecod=%d ener=%f",
                 idcode[i], icou[i], imod[i], z[i], a[i], ecode[i], ener[i]);
         continue;
      }
      auto node = det->GetNode();
      auto traj = (KVGeoDNTrajectory*)node->GetTrajectories()->First();
      bool etalon = sili_fired(i) || si75_fired(i);
      auto chio = det->GetChIo();
      if (node->GetNTrajForwards() > 1) {
         // Stopped in CsI behind etalon telescope - we have a particle stopped in CsI which may have gone CI-CsI or CI-Si75-SiLi-CsI
         if (etalon) {
            if (stopped_in_csi(i)) {
               // particles stopping in CsI can only be on etalon trajectory if they have energy necessary
               // to punch through CI-SI75-SILI detectors
               KVDetectorStack etalon_dets;
               etalon_dets.Add(chio);
               etalon_dets.Add(get_si75(i));
               etalon_dets.Add(get_sili(i));
               if (z[i] > 0 && z[i] < 92 && std::abs(a[i]) > 0 && ecode[i] > 0 && ecode[i] < 4 && ener[i] > 0) {
                  if (ener[i] < etalon_dets.GetPunchThroughEnergy(z[i], std::abs(a[i]))) {
                     etalon = false;
                     traj = (KVGeoDNTrajectory*)node->GetTrajectories()->FindObjectWithMethod("2", "GetN");
                  }
               }
               else
                  traj = (KVGeoDNTrajectory*)node->GetTrajectories()->FindObjectWithMethod("4", "GetN");
            }
            else
               traj = (KVGeoDNTrajectory*)node->GetTrajectories()->FindObjectWithMethod("4", "GetN");
         }
         else
            traj = (KVGeoDNTrajectory*)node->GetTrajectories()->FindObjectWithMethod("2", "GetN");
      }
      // correct idtelescope code in case of code=6 or code=10 stopped in SiLi or Si75 (1st campaign data)
      if (etalon && (idcode[i] == KVINDRA::ID_CI_SI_COHERENCY || idcode[i] == KVINDRA::ID_CSI_MASS_OUT_OF_RANGE)) {
         if (stopped_in_sili(i))
            idtelescope_code = KVINDRA::ID_SI75_SILI;
         else if (stopped_in_si75(i))
            idtelescope_code = KVINDRA::ID_CI_SI75;
      }
 
      if (!node->GetTrajectories()) {
         Warning("Process", "node %s has no trajectories", node->GetName());
      }
      auto rtraj = (KVReconNucTrajectory*)det->GetGroup()->GetTrajectoryForReconstruction(traj, node);
      TIter next_id(rtraj->GetIDTelescopes());
      KVIDTelescope* idt = nullptr;
      while ((idt = (KVIDTelescope*)next_id())) {
         if (idt->GetIDCode() == idtelescope_code) break;
      }
      if (!idt || idt->GetIDCode() != idtelescope_code) {
         idt = (KVIDTelescope*)rtraj->GetIDTelescopes()->First();
      }
      // idt may still = nullptr here if e.g. particle stopped in ChIo (idcode=stopped_in_first_stage)
      // i.e. normally only in cases where no identification is available
 
      auto nuc = the_event->AddParticle();
      nuc->SetParameter("ARRAY", "INDRA");
      nuc->SetMassFormula(KVNucleus::kVedaMass);
      nuc->SetReconstructionTrajectory(rtraj);
      if (phos_fired(i)) {
         the_event->SetParameter(Form("ACQPAR.INDRA.PHOS_%02d_R", imod[i]), ci_gg[i]);
         the_event->SetParameter(Form("ACQPAR.INDRA.PHOS_%02d_L", imod[i]), ci_pg[i]);
      }
      if (ci_fired(i)) {
         if (chio) {
            the_event->SetParameter(Form("ACQPAR.INDRA.%s_GG", chio->GetName()), ci_gg[i]);
            the_event->SetParameter(Form("ACQPAR.INDRA.%s_PG", chio->GetName()), ci_pg[i]);
         }
      };
      if (si_fired(i)) {
         the_event->SetParameter(Form("ACQPAR.INDRA.SI_%02d%02d_GG", icou[i], imod[i]), si_gg[i]);
         the_event->SetParameter(Form("ACQPAR.INDRA.SI_%02d%02d_PG", icou[i], imod[i]), si_pg[i]);
      }
      if (csi_fired(i)) {
         the_event->SetParameter(Form("ACQPAR.INDRA.CSI_%02d%02d_R", icou[i], imod[i]), csi_r[i]);
         the_event->SetParameter(Form("ACQPAR.INDRA.CSI_%02d%02d_L", icou[i], imod[i]), csi_l[i]);
      }
 
      // set up identification result
      auto idr = nuc->GetIdentificationResult(1);
      idr->IDattempted = true;
      if (idt) {
         idr->SetIDType(idt->GetType());
         idr->IDOK = true;
         idr->Zident = true;
      }
      idr->IDcode = idcode[i];
      // correction for 1st campaign data: ALL masses are positive!!!
      if (camp1) {
         // after inspection of data, only CsI R/L provides isotopic identification (for Z<5)
         if (idcode[i] != 2) {
            a_indra[i] = -TMath::Abs(a_indra[i]);
            a[i] = -TMath::Abs(a[i]);
         }
      }
      idr->Aident = (a[i] > 0);
      if (idr->Aident) idr->PID = a_indra[i];
      else idr->PID = z_indra[i];
      if (z[i] < 0 || z[i] > 92)
         idr->Z = 0;
      else
         idr->Z = z[i];
      idr->A = TMath::Abs(a[i]);
      if (idcode[i] != KVINDRA::ID_CI_COHERENCY
            && idcode[i] != KVINDRA::ID_CI_MULTIHIT
            && idcode[i] != KVINDRA::ID_CSI_FRAGMENT
            && idcode[i] != KVINDRA::ID_CSI_MASS_OUT_OF_RANGE
            && idcode[i] != KVINDRA::ID_STOPPED_IN_FIRST_STAGE)
         nuc->SetIsIdentified();
      nuc->SetIdentification(idr, idt);
 
      if (ecode[i] > 0 && ecode[i] < 4 && ener[i] > 0) {
         nuc->SetIsCalibrated();
         nuc->SetECode(ecode[i]);
         nuc->SetEnergy(ener[i]);
         nuc->GetAnglesFromReconstructionTrajectory();
         if (phos_fired(i)) {
            nuc->SetParameter("INDRA.EPHOS", de1[i]);
         }
         if (ci_fired(i)) {
            nuc->SetParameter("INDRA.ECHIO", de1[i]);
            nuc->SetParameter("INDRA.DE_MYLAR", de_mylar[i]);
         };
         if (si_fired(i)) {
            nuc->SetParameter("INDRA.ESI", de2[i]);
         }
         if (csi_fired(i)) {
            nuc->SetParameter("INDRA.ECSI", de3[i]);
         }
         if (si75_fired(i)) nuc->SetParameter("INDRA.ESI75", de4[i]);
         if (sili_fired(i)) nuc->SetParameter("INDRA.ESILI", de5[i]);
      }
      else
         nuc->SetIsUncalibrated();
   }
 
   output_tree->Fill();
 
   ++events_read;
   if (!(events_read % 1000)) std::cout << "Read " << events_read << "/" << total_events << " events...\n";
 
   return kTRUE;
}
 
 
 
 
void KVINDRADSTConverter::SlaveTerminate()
{
   // The SlaveTerminate() function is called after all entries or objects
   // have been processed. When running with PROOF SlaveTerminate() is called
   // on each slave server.
 
}
 
 
 
 
void KVINDRADSTConverter::Terminate()
{
   // The Terminate() function is the last function to be called during
   // a query. It always runs on the client, it can be used to present
   // the results graphically or save the results to file.
 
   output_file->cd();
   gDataAnalyser->WriteBatchInfo(output_tree);
   output_tree->Write();     //write tree to file
 
   //add new file to repository
   gDataSet->CommitRunfile("recon", run_index, output_file);
}