diff --git a/macro/littrack/radlength_ana.C b/macro/littrack/radlength_ana.C
index 314dae92697a4e9132f64426c92577dcba5b61e3..4957c7581bab35a647b30b2a25efde21697169c7 100644
--- a/macro/littrack/radlength_ana.C
+++ b/macro/littrack/radlength_ana.C
@@ -4,58 +4,69 @@
* \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"));
- }
-
- TStopwatch timer;
- timer.Start();
-
- gROOT->LoadMacro("$VMCWORKDIR/macro/littrack/loadlibs.C");
- loadlibs();
-
- // ----- Reconstruction run -------------------------------------------
- FairRunAna* run = new FairRunAna();
- run->SetInputFile(mcFile);
- run->SetOutputFile(radLengthQaFile);
- // ------------------------------------------------------------------------
-
- CbmLitRadLengthQa* radLengthQa = new CbmLitRadLengthQa();
- radLengthQa->SetOutputDir(std::string(resultDir));
- run->AddTask(radLengthQa);
-
- // ----- Parameter database --------------------------------------------
- FairRuntimeDb* rtdb = run->GetRuntimeDb();
- FairParRootFileIo* parIo1 = new FairParRootFileIo();
- parIo1->open(parFile.Data());
- rtdb->setFirstInput(parIo1);
- rtdb->setOutput(parIo1);
- rtdb->saveOutput();
- // ------------------------------------------------------------------------
-
- // ----- 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 << "Parameter file is " << parFile << std::endl;
- std::cout << "Real time " << timer.RealTime() << " s, CPU time "
- << timer.CpuTime() << " s" << std::endl;
- // ------------------------------------------------------------------------
+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;
+
+ FairRunAna *run = new FairRunAna();
+ FairFileSource* inputSource = new FairFileSource(mcFile.c_str());
+ run->SetSource(inputSource);
+ run->SetOutputFile(radqaFile.c_str());
+ run->SetGenerateRunInfo(kTRUE);
+
+
+ CbmLitRadLengthQa* radLengthQa = new CbmLitRadLengthQa();
+ radLengthQa->SetOutputDir(resultDir);
+ run->AddTask(radLengthQa);
+
+ 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();
+
+ run->Run(0, nEvents);
+
+ timer.Stop();
+ 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;
}
diff --git a/macro/littrack/radlength_sim.C b/macro/littrack/radlength_sim.C
index d22ab81140a5d2b8dc4936a6a8432ea7c35c8d41..316fc1b31e16ed89a9ba4ad1248c925ae3c4e065 100644
--- a/macro/littrack/radlength_sim.C
+++ b/macro/littrack/radlength_sim.C
@@ -3,99 +3,37 @@
using std::cout;
using std::endl;
-void radlength_sim(Int_t nEvents = 1210000) {
- TString script = TString(gSystem->Getenv("LIT_SCRIPT"));
+void radlength_sim(
+ 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 &geoSetup = "sis100_electron_rich_pal_bcarb",
+ Int_t nofEvents = 0) {
- TString dir = "events/radlength_rich_v13c/";
- TString mcFile = dir + "radlength.mc.0000.root";
- TString parFile = dir + "radlength.param.0000.root";
+ TTree::SetMaxTreeSize(90000000000);
- TString caveGeom = "cave.geo";
- TString mvdGeom = ""; //"mvd/mvd_v07a.geo";
- TString stsGeom = ""; //"sts/sts_v12b.geo.root";
- TString richGeom = "rich/rich_v13c.root";
- TString trdGeom = ""; //"trd/trd_v13p_3e.geo.root";
- TString muchGeom = ""; //"much/much_v12b.geo";
- TString tofGeom = ""; //"tof/tof_v13b.geo.root";
-
- if (script == "yes") {
- mcFile = TString(gSystem->Getenv("LIT_MC_FILE"));
- parFile = TString(gSystem->Getenv("LIT_PAR_FILE"));
-
- caveGeom = TString(gSystem->Getenv("LIT_CAVE_GEOM"));
- //stsGeom = TString(gSystem->Getenv("LIT_STS_GEOM"));
- trdGeom = TString(gSystem->Getenv("LIT_TRD_GEOM"));
- }
+ remove(parFile.c_str());
+ remove(mcFile.c_str());
TStopwatch timer;
timer.Start();
- gROOT->LoadMacro("$VMCWORKDIR/macro/littrack/loadlibs.C");
- loadlibs();
-
- FairRunSim* run = new FairRunSim();
- run->SetName("TGeant3"); // Transport engine
- run->SetOutputFile(mcFile); // Output file
- FairRuntimeDb* rtdb = run->GetRuntimeDb();
-
- run->SetMaterials("media.geo"); // Materials
-
- // ----- Create detectors and passive volumes -------------------------
- if (caveGeom != "") {
- FairModule* cave = new CbmCave("CAVE");
- cave->SetGeometryFileName(caveGeom);
- run->AddModule(cave);
- }
-
- if (mvdGeom != "") {
- FairDetector* mvd = new CbmMvd("MVD", kTRUE);
- mvd->SetGeometryFileName(mvdGeom);
- run->AddModule(mvd);
- }
-
- if (stsGeom != "") {
- FairDetector* sts = new CbmStsMC(kTRUE);
- sts->SetGeometryFileName(stsGeom);
- run->AddModule(sts);
- }
-
- if (richGeom != "") {
- FairDetector* rich = new CbmRich("RICH", kTRUE);
- rich->SetGeometryFileName(richGeom);
- run->AddModule(rich);
- }
-
- if (muchGeom != "") {
- FairDetector* much = new CbmMuch("MUCH", kTRUE);
- much->SetGeometryFileName(muchGeom);
- run->AddModule(much);
- }
-
- if (trdGeom != "") {
- FairDetector* trd = new CbmTrd("TRD", kTRUE);
- trd->SetGeometryFileName(trdGeom);
- run->AddModule(trd);
- }
-
- if (tofGeom != "") {
- FairDetector* tof = new CbmTof("TOF", kTRUE);
- tof->SetGeometryFileName(tofGeom);
- run->AddModule(tof);
- }
-
- // ----- Create PrimaryGenerator --------------------------------------
- FairPrimaryGenerator* primGen = new FairPrimaryGenerator();
- run->SetGenerator(primGen);
+// const Double_t minX = -550; // cm
+// const Double_t maxX = 550; // cm
+// const Double_t stepX = 10.; // cm
+// const Int_t nofBinsX = (maxX - minX) / stepX;
+// const Double_t minY = -550; // cm
+// const Double_t maxY = 550; // cm
+// const Double_t stepY = 10.; // cm
- const Double_t minX = -550; // cm
- const Double_t maxX = 550; // cm
+ const Double_t minX = -250; // cm
+ const Double_t maxX = 250; // cm
const Double_t stepX = 1.; // cm
const Int_t nofBinsX = (maxX - minX) / stepX;
- const Double_t minY = -550; // cm
- const Double_t maxY = 550; // cm
+ const Double_t minY = -250; // cm
+ const Double_t maxY = 250; // cm
const Double_t stepY = 1.; // cm
const Int_t nofBinsY = (maxY - minY) / stepY;
- const Int_t nofEvents = nofBinsX * nofBinsY;
+ nofEvents = nofBinsX * nofBinsY;
std::vector<Double_t> vectorX, vectorY;
for (Int_t iX = 0; iX < nofBinsX; iX++) {
@@ -107,13 +45,10 @@ void radlength_sim(Int_t nEvents = 1210000) {
}
}
- std::cout << "CbmLitRadLengthGenerator:\n";
- std::cout << "nofBinX=" << nofBinsX << " nofBinsY=" << nofBinsY
- << " nofEvents=" << nofEvents << std::endl;
-
+ CbmTransport run;
CbmLitRadLengthGenerator* generator = new CbmLitRadLengthGenerator();
generator->SetXY(vectorX, vectorY);
- primGen->AddGenerator(generator);
+ run.AddInput(generator);
/*
const int RMax = 700; // Maximum radius of the station
@@ -133,29 +68,22 @@ void radlength_sim(Int_t nEvents = 1210000) {
box->SetThetaRange(0., 50.);
primGen->AddGenerator(box);
*/
- run->SetStoreTraj(kFALSE);
- run->SetRadLenRegister(kTRUE);
- // ------------------------------------------------------------------------
-
- run->Init();
-
- // ----- Runtime database ---------------------------------------------
- Bool_t kParameterMerged = kTRUE;
- FairParRootFileIo* parOut = new FairParRootFileIo(kParameterMerged);
- parOut->open(parFile.Data());
- rtdb->setOutput(parOut);
- rtdb->saveOutput();
- rtdb->print();
- // ------------------------------------------------------------------------
-
- run->Run(nEvents);
- // ----- Finish -------------------------------------------------------
- timer.Stop();
- cout << "Macro finished succesfully." << endl;
- cout << "Output file is " << mcFile << endl;
- cout << "Parameter file is " << parFile << endl;
- cout << "Real time " << timer.RealTime() << " s, CPU time " << timer.CpuTime()
- << "s" << endl;
- // ------------------------------------------------------------------------
+ run.SetOutFileName(mcFile.c_str());
+ run.SetParFileName(parFile.c_str());
+ run.LoadSetup(geoSetup.c_str());
+ run.SetTarget("Gold", 0.025, 2.5);
+ run.SetBeamPosition(0., 0., 0., 0.);
+ //run.SetEngine(kGeant4);
+ //run.StoreTrajectories(true);
+ run.RegisterRadLength(true);
+ run.Run(nofEvents);
+
+ timer.Stop();
+ std::cout << std::endl << std::endl;
+ std::cout << "Macro finished successfully." << std::endl;
+ std::cout << "MC file is " << mcFile << std::endl;
+ std::cout << "Parameter file is " << parFile << std::endl;
+ std::cout << "Real time " << timer.RealTime() << " s, CPU time " << timer.CpuTime() << "s" << std::endl;
+ std::cout << std::endl << "Test passed" << std::endl << "All ok" << std::endl;
}
diff --git a/reco/littrack/cbm/qa/radlength/CbmLitRadLengthQa.cxx b/reco/littrack/cbm/qa/radlength/CbmLitRadLengthQa.cxx
index 059d1f71d77817b3fca9fb8ce2caae8525d93fe8..3637541572aa7e877d5b882ba87b39b3850e3bfb 100644
--- a/reco/littrack/cbm/qa/radlength/CbmLitRadLengthQa.cxx
+++ b/reco/littrack/cbm/qa/radlength/CbmLitRadLengthQa.cxx
@@ -63,17 +63,17 @@ void CbmLitRadLengthQa::Exec(Option_t* opt) {
}
ExecDetector(".+", "Total");
- ExecDetector("/cave_1/pipevac1_0/mvdstation.+", "Mvd");
- ExecDetector("/cave_1/STS.+", "Sts");
+ //ExecDetector("/cave_1/pipevac1_0/mvdstation.+", "Mvd");
+ ExecDetector("/cave_1/sts.+", "Sts");
ExecDetector("/cave_1/rich.+", "Rich");
- ExecDetector("/cave_1/trd.+", "Trd");
- ExecDetector("/cave_1/much.+", "Much");
- ExecDetector("/cave_1/tof.+", "Tof");
- ExecDetector("Mvd", CbmLitRadLengthQa::GetMvdStationId);
- ExecDetector("Sts", CbmLitRadLengthQa::GetStsStationId);
- ExecDetector("Trd", CbmLitRadLengthQa::GetTrdStationId);
- ExecDetector("Much", CbmLitRadLengthQa::GetMuchStationId);
- ExecDetector("MuchAbsorber", CbmLitRadLengthQa::GetMuchAbsorberId);
+ //ExecDetector("/cave_1/trd.+", "Trd");
+ //ExecDetector("/cave_1/much.+", "Much");
+ //ExecDetector("/cave_1/tof.+", "Tof");
+ //ExecDetector("Mvd", CbmLitRadLengthQa::GetMvdStationId);
+ //ExecDetector("Sts", CbmLitRadLengthQa::GetStsStationId);
+ //ExecDetector("Trd", CbmLitRadLengthQa::GetTrdStationId);
+ //ExecDetector("Much", CbmLitRadLengthQa::GetMuchStationId);
+ //ExecDetector("MuchAbsorber", CbmLitRadLengthQa::GetMuchAbsorberId);
}
void CbmLitRadLengthQa::Finish() {
@@ -100,15 +100,15 @@ void CbmLitRadLengthQa::ReadDataBranches() {
}
void CbmLitRadLengthQa::CreateHistograms() {
- const Int_t nofBins = 1000;
- const Int_t nofBinsX = 1000;
- const Int_t nofBinsY = 1000;
- const Int_t nofBinsSiliconThicknessX = 1100;
- const Int_t nofBinsSiliconThicknessY = 1100;
- const Double_t minX = -550.1;
- const Double_t maxX = 549.9;
- const Double_t minY = -550.1;
- const Double_t maxY = 549.9;
+ const Int_t nofBins = 500;
+ const Int_t nofBinsX = 500;
+ const Int_t nofBinsY = 500;
+ const Int_t nofBinsSiliconThicknessX = 500;
+ const Int_t nofBinsSiliconThicknessY = 500;
+ const Double_t minX = -250.1;
+ const Double_t maxX = 249.9;
+ const Double_t minY = -250.1;
+ const Double_t maxY = 249.9;
vector<string> detNames =
list_of("Total")("Mvd")("Sts")("Rich")("Trd")("Much")("Tof");
Int_t nofDetNames = detNames.size();
diff --git a/reco/littrack/cbm/qa/radlength/CbmLitRadLengthQaReport.cxx b/reco/littrack/cbm/qa/radlength/CbmLitRadLengthQaReport.cxx
index 56f6291afd9d0cecdcf35c9e0f6de5c767bf6f77..fb7f1c583002b08a8e6e019112fa74678da44fc1 100644
--- a/reco/littrack/cbm/qa/radlength/CbmLitRadLengthQaReport.cxx
+++ b/reco/littrack/cbm/qa/radlength/CbmLitRadLengthQaReport.cxx
@@ -27,6 +27,8 @@ void CbmLitRadLengthQaReport::Create() {
}
void CbmLitRadLengthQaReport::Draw() {
+ SetDefaultDrawStyle();
+
DrawDetector("Total");
DrawDetector("Mvd");
DrawDetector("Sts");
@@ -77,6 +79,17 @@ void CbmLitRadLengthQaReport::DrawDetector(const string& detName) {
1200,
1000);
DrawH2(HM()->P2("hrl_RadThickness_" + detName + "_P2"));
+
+
+ TCanvas* canvas4_2 =
+ CreateCanvas(string("hrl_RadThickness_" + detName + "_P2_zoom").c_str(),
+ string("hrl_RadThickness_" + detName + "_P2_zoom").c_str(),
+ 1200,
+ 1000);
+ TProfile2D* pr =
+ (TProfile2D*) HM()->P2("hrl_RadThickness_" + detName + "_P2")->Clone();
+ pr->GetZaxis()->SetRangeUser(0, 100);
+ DrawH2(pr);
}
if (HM()->Exists("hrl_ThicknessSilicon_" + detName + "_H1")) {