Updates in radiation length calculation

macro/littrack/radlength_ana.C

@@ -4,58 +4,72 @@
  * \author Andrey Lebedev <andrey.lebedev@gsi.de>
  * \date 2013
  */
-void radlength_ana(Int_t nEvents = 1210000) {
-  TString script = TString(gSystem->Getenv("LIT_SCRIPT"));
-
-  TString dir     = "events/radlength_rich_v13c/";      // Output directory
-  TString mcFile  = dir + "radlength.mc.0000.root";     // MC transport file
-  TString parFile = dir + "radlength.param.0000.root";  // Parameter file
-  TString radLengthQaFile = dir + "radlength.qa.root";  // MC transport file
-  TString resultDir       = "events/results_radlength_rich_v13c/";
-
-  if (script == "yes") {
-    mcFile          = TString(gSystem->Getenv("LIT_MC_FILE"));
-    parFile         = TString(gSystem->Getenv("LIT_PAR_FILE"));
-    radLengthQaFile = TString(gSystem->Getenv("LIT_RAD_LENGTH_QA_FILE"));
-    resultDir       = TString(gSystem->Getenv("LIT_RESULT_DIR"));
-  }
+void radlength_ana(
+  const string& mcFile =
+    "/Users/slebedev/Development/cbm/data/sim/rich/radlen/mc.ac.root",
+  const string& parFile =
+    "/Users/slebedev/Development/cbm/data/sim/rich/radlen/param.ac.root",
+  const string& radqaFile =
+    "/Users/slebedev/Development/cbm/data/sim/rich/radlen/radqa.ac.root",
+  const string& resultDir = "results_radlen_ac/",
+  const string& geoSetup  = "sis100_electron_rich_pal_bcarb",
+  Int_t nEvents           = 12100000) {
+
+  FairLogger::GetLogger()->SetLogScreenLevel("INFO");
+  FairLogger::GetLogger()->SetLogVerbosityLevel("LOW");
+  TTree::SetMaxTreeSize(90000000000);
+
+  TString srcDir = gSystem->Getenv("VMCWORKDIR");  // top source directory
+
+  remove(radqaFile.c_str());
+
+  TString setupFile  = srcDir + "/geometry/setup/setup_" + geoSetup + ".C";
+  TString setupFunct = "setup_" + geoSetup + "()";
+  gROOT->LoadMacro(setupFile);
+  gROOT->ProcessLine(setupFunct);
+
+  TList* parFileList = new TList();
 
   TStopwatch timer;
   timer.Start();
+  gDebug = 0;
 
-  gROOT->LoadMacro("$VMCWORKDIR/macro/littrack/loadlibs.C");
-  loadlibs();
+  FairRunAna* run             = new FairRunAna();
+  FairFileSource* inputSource = new FairFileSource(mcFile.c_str());
+  run->SetSource(inputSource);
+  run->SetOutputFile(radqaFile.c_str());
+  run->SetGenerateRunInfo(kTRUE);
 
-  // -----   Reconstruction run   -------------------------------------------
-  FairRunAna* run = new FairRunAna();
-  run->SetInputFile(mcFile);
-  run->SetOutputFile(radLengthQaFile);
-  // ------------------------------------------------------------------------
 
   CbmLitRadLengthQa* radLengthQa = new CbmLitRadLengthQa();
-  radLengthQa->SetOutputDir(std::string(resultDir));
+  radLengthQa->SetOutputDir(resultDir);
   run->AddTask(radLengthQa);
 
-  // -----  Parameter database   --------------------------------------------
-  FairRuntimeDb* rtdb       = run->GetRuntimeDb();
-  FairParRootFileIo* parIo1 = new FairParRootFileIo();
-  parIo1->open(parFile.Data());
+  FairRuntimeDb* rtdb        = run->GetRuntimeDb();
+  FairParRootFileIo* parIo1  = new FairParRootFileIo();
+  FairParAsciiFileIo* parIo2 = new FairParAsciiFileIo();
+  parIo1->open(parFile.c_str(), "UPDATE");
   rtdb->setFirstInput(parIo1);
+  if (!parFileList->IsEmpty()) {
+    parIo2->open(parFileList, "in");
+    rtdb->setSecondInput(parIo2);
+  }
+
+  run->Init();
+
   rtdb->setOutput(parIo1);
   rtdb->saveOutput();
-  // ------------------------------------------------------------------------
+  rtdb->print();
 
-  // -----   Initialize and run   --------------------------------------------
-  run->Init();
   run->Run(0, nEvents);
-  // ------------------------------------------------------------------------
 
-  // -----   Finish   -------------------------------------------------------
   timer.Stop();
-  std::cout << "Macro finished successfully." << std::endl;
-  std::cout << "Output file is " << radLengthQaFile << std::endl;
+  std::cout << std::endl << std::endl;
+  std::cout << "Macro finished succesfully." << std::endl;
+  std::cout << "Output file is " << mcFile << std::endl;
   std::cout << "Parameter file is " << parFile << std::endl;
+  std::cout << "Radqa file is " << radqaFile << std::endl;
   std::cout << "Real time " << timer.RealTime() << " s, CPU time "
             << timer.CpuTime() << " s" << std::endl;
-  // ------------------------------------------------------------------------
+  std::cout << "Test passed" << std::endl << "All ok" << std::endl;
 }