diff --git a/macro/alignment/run_BbaAlignment_mcbm.C b/macro/alignment/run_BbaAlignment_mcbm.C
index 66dbc65a52c074a2e03bc2ef490018c0b101c7c6..dedcce7ae56e3f905a3be5ca02cebb8fcf7a13e0 100644
--- a/macro/alignment/run_BbaAlignment_mcbm.C
+++ b/macro/alignment/run_BbaAlignment_mcbm.C
@@ -36,8 +36,9 @@
#include <TStopwatch.h>
#endif
-void run_BbaAlignment_mcbm(Int_t nEvents = 0, TString dataset = "data/mcbm_beam_2020_03_test",
- TString setupName = "mcbm_beam_2020_03", Bool_t bUseMC = kTRUE, const TString& config = "")
+void run_BbaAlignment_mcbm(Int_t nEvents = -1, TString dataset = "data/mcbm_beam_2020_03_test",
+ TString setupName = "mcbm_beam_2022_05_23_nickel", double SimulatedMisalignmentRange = 0.,
+ Bool_t bUseMC = kFALSE)
{
// ========================================================================
@@ -49,92 +50,27 @@ void run_BbaAlignment_mcbm(Int_t nEvents = 0, TString dataset = "data/mcbm_beam_
// ------------------------------------------------------------------------
// ----- Environment --------------------------------------------------
- int verbose = 6; // verbose level
- TString myName = "mcbm_qa"; // this macro's name for screen output
- TString srcDir = gSystem->Getenv("VMCWORKDIR"); // top source directory
- TString qaConfig = (config.Length() ? config : srcDir + "/macro/qa/configs/qa_tasks_config_mcbm.yaml");
+ int verbose = 6; // verbose level
+ TString myName = "BBA Alignment"; // this macro's name for screen output
+ TString srcDir = gSystem->Getenv("VMCWORKDIR"); // top source directory
// NOTE: SZh 28.07.2024: config can depend from the setup
// ------------------------------------------------------------------------
// ----- In- and output file names ------------------------------------
- TString rawFile = dataset + ".event.raw.root";
- TString traFile = dataset + ".tra.root";
+ TString traFile = dataset + ".tra.root";
+ //TString rawFile = dataset + ".event.raw.root";
TString parFile = dataset + ".par.root";
TString recFile = dataset + ".rec.root";
- TString sinkFile = dataset + ".qa.root";
+ TString sinkFile = dataset + ".ali.root";
+ TString geoFile = dataset + ".geo.root";
+
// ------------------------------------------------------------------------
// ----- Load the geometry setup -------------------------------------
std::cout << '\n';
- TString setupFile = srcDir + "/geometry/setup/setup_" + setupName + ".C";
- TString setupFunct = "setup_";
- setupFunct = setupFunct + setupName + "()";
- std::cout << "-I- " << myName << ": Loading macro " << setupFile << '\n';
- gROOT->LoadMacro(setupFile);
- gROOT->ProcessLine(setupFunct);
- CbmSetup* setup = CbmSetup::Instance();
- // setup->RemoveModule(ECbmModuleId::kTrd);
- // ------------------------------------------------------------------------
- // ----- Some global switches -----------------------------------------
- //if (setupName == "mcbm_beam_2022_05_23_nickel") { setup->RemoveModule(ECbmModuleId::kMuch); }
-
- //bool eventBased = !sEvBuildRaw.IsNull();
- bool bUseMvd = setup->IsActive(ECbmModuleId::kMvd);
- bool bUseSts = setup->IsActive(ECbmModuleId::kSts);
- bool bUseRich = setup->IsActive(ECbmModuleId::kRich);
- bool bUseMuch = setup->IsActive(ECbmModuleId::kMuch);
- bool bUseTrd = setup->IsActive(ECbmModuleId::kTrd);
- bool bUseTof = setup->IsActive(ECbmModuleId::kTof);
- bool bUsePsd = setup->IsActive(ECbmModuleId::kPsd);
- std::cout << " MVD: " << (bUseMvd ? "ON" : "OFF") << '\n';
- std::cout << " STS: " << (bUseSts ? "ON" : "OFF") << '\n';
- std::cout << " RICH: " << (bUseRich ? "ON" : "OFF") << '\n';
- std::cout << " MUCH: " << (bUseMuch ? "ON" : "OFF") << '\n';
- std::cout << " TRD: " << (bUseTrd ? "ON" : "OFF") << '\n';
- std::cout << " TOF: " << (bUseTof ? "ON" : "OFF") << '\n';
- std::cout << " PSD: " << (bUsePsd ? "ON" : "OFF") << '\n';
- // ------------------------------------------------------------------------
-
- // ----- Parameter files as input to the runtime database -------------
- std::cout << '\n';
- std::cout << "-I- " << myName << ": Defining paramete files\n";
- TList* parFileList = new TList();
- TString geoTag;
-
- // - MUCH digitisation parameters
- TString muchParFile{};
- if (setup->GetGeoTag(ECbmModuleId::kMuch, geoTag)) {
- bool mcbmFlag = geoTag.Contains("mcbm", TString::kIgnoreCase);
- muchParFile = srcDir + "/parameters/much/much_";
- muchParFile += (mcbmFlag) ? geoTag : geoTag(0, 4);
- muchParFile += "_digi_sector.root";
- { // init geometry from the file
- TFile* f = new TFile(muchParFile, "R");
- TObjArray* stations = f->Get<TObjArray>("stations");
- assert(stations);
- CbmMuchGeoScheme::Instance()->Init(stations, mcbmFlag);
- }
- }
-
- // - TRD digitisation parameters
- if (setup->GetGeoTag(ECbmModuleId::kTrd, geoTag)) {
- const Char_t* npar[4] = {"asic", "digi", "gas", "gain"};
- TObjString* trdParFile(NULL);
- for (Int_t i(0); i < 4; i++) {
- trdParFile = new TObjString(srcDir + "/parameters/trd/trd_" + geoTag + "." + npar[i] + ".par");
- parFileList->Add(trdParFile);
- std::cout << "-I- " << myName << ": Using parameter file " << trdParFile->GetString() << std::endl;
- }
- }
-
- // - TOF digitisation parameters
- if (setup->GetGeoTag(ECbmModuleId::kTof, geoTag)) {
- TObjString* tofBdfFile = new TObjString(srcDir + "/parameters/tof/tof_" + geoTag + ".digibdf.par");
- parFileList->Add(tofBdfFile);
- std::cout << "-I- " << myName << ": Using parameter file " << tofBdfFile->GetString() << std::endl;
- }
- // ------------------------------------------------------------------------
+ CbmSetup* setup = CbmSetup::Instance();
+ setup->LoadSetup(setupName);
// In general, the following parts need not be touched
// ========================================================================
@@ -149,14 +85,16 @@ void run_BbaAlignment_mcbm(Int_t nEvents = 0, TString dataset = "data/mcbm_beam_
// ------------------------------------------------------------------------
// ----- FairRunAna ---------------------------------------------------
- FairFileSource* inputSource = new FairFileSource(rawFile);
- inputSource->AddFriend(traFile);
- inputSource->AddFriend(recFile);
+ FairFileSource* inputSource = new FairFileSource(recFile);
+ if (bUseMC) {
+ inputSource->AddFriend(traFile);
+ }
+ //inputSource->AddFriend(rawFile);
FairRunAna* run = new FairRunAna();
run->SetSource(inputSource);
run->SetGenerateRunInfo(kFALSE);
-
+ run->SetGeomFile(geoFile);
FairRootFileSink* sink = new FairRootFileSink(sinkFile);
run->SetSink(sink);
@@ -187,11 +125,6 @@ void run_BbaAlignment_mcbm(Int_t nEvents = 0, TString dataset = "data/mcbm_beam_
// ------------------------------------------------------------------------
- // ----- CA tracking QA ---------------------------------------------------
- // Kalman Filter (currently needed to access the magnetic filed, to be
- // removed soon)
- run->AddTask(new CbmKF());
-
// ------------------------------------------------------------------------
// TODO: read tracking parameters from the file like CaOutputQa does
@@ -212,8 +145,10 @@ void run_BbaAlignment_mcbm(Int_t nEvents = 0, TString dataset = "data/mcbm_beam_
// ----- BBA alignment --------------------------------------------
CbmBbaAlignmentTask* alignment = new CbmBbaAlignmentTask();
alignment->SetMcbmTrackingMode();
- run->AddTask(alignment);
+ alignment->SetSimulatedMisalignmentRange(SimulatedMisalignmentRange);
+ run->AddTask(alignment);
+ run->SetGenerateRunInfo(kFALSE);
// ----- Parameter database --------------------------------------------
std::cout << std::endl << std::endl;
@@ -222,12 +157,7 @@ void run_BbaAlignment_mcbm(Int_t nEvents = 0, TString dataset = "data/mcbm_beam_
FairParRootFileIo* parIo1 = new FairParRootFileIo();
parIo1->open(parFile.Data(), "in");
rtdb->setFirstInput(parIo1);
- if (!parFileList->IsEmpty()) {
- FairParAsciiFileIo* parIo2 = new FairParAsciiFileIo();
- parIo2->open(parFileList, "in");
- rtdb->setSecondInput(parIo2);
- }
- rtdb->print();
+
// ------------------------------------------------------------------------
// ----- Run initialisation -------------------------------------------
@@ -235,6 +165,8 @@ void run_BbaAlignment_mcbm(Int_t nEvents = 0, TString dataset = "data/mcbm_beam_
std::cout << "-I- " << myName << ": Initialise run" << std::endl;
run->Init();
+ rtdb->print();
+
// ----- Start run ----------------------------------------------------
std::cout << std::endl << std::endl;
std::cout << "-I- " << myName << ": Starting run" << std::endl;
@@ -244,6 +176,7 @@ void run_BbaAlignment_mcbm(Int_t nEvents = 0, TString dataset = "data/mcbm_beam_
// ----- Finish -------------------------------------------------------
timer.Stop();
+
FairMonitor::GetMonitor()->Print();
Double_t rtime = timer.RealTime();
Double_t ctime = timer.CpuTime();
diff --git a/macro/beamtime/mcbm2022/mcbm_event_reco_L1.C b/macro/beamtime/mcbm2022/mcbm_event_reco_L1.C
index 6494eb9a84fc41837399072376fd76b334493b67..9fdb67ee0c01de3d6743cb23c98256cbecf8d357 100644
--- a/macro/beamtime/mcbm2022/mcbm_event_reco_L1.C
+++ b/macro/beamtime/mcbm2022/mcbm_event_reco_L1.C
@@ -59,9 +59,11 @@ Bool_t mcbm_event_reco_L1(UInt_t uRunId = 2570,
TString inFile = sInpDir + sRunId + ".digi";
TString cFileId = sRunId + ".5.000";
+ gSystem->Exec("mkdir " + sOutDir);
+ gSystem->Exec("cp $VMCWORKDIR/macro/run/.rootrc .");
+
//TString parFileIn = sInpDir + "/unp_mcbm_params_" + sRunId;
- TString parFileOut = sOutDir + "/reco_event_mcbm_test_params_" + sRunId;
- TString outFile = sOutDir + "/reco_event_mcbm_test" + sRunId;
+ TString outFile = sOutDir + "/mcbm_test" + sRunId;
// Your folder with the Tof Calibration files;
TString TofFileFolder = "";
@@ -72,14 +74,16 @@ Bool_t mcbm_event_reco_L1(UInt_t uRunId = 2570,
if (0 <= iUnpFileIndex) {
inFile += TString::Format("_%02u", iUnpFileIndex);
// the input par file is not split during unpacking!
- parFileOut += TString::Format("_%02u", iUnpFileIndex);
outFile += TString::Format("_%02u", iUnpFileIndex);
} // if ( 0 <= uUnpFileIndex )
/// Add ROOT file suffix
inFile += ".root";
// parFileIn += ".root";
- parFileOut += ".root";
- outFile += ".root";
+
+ TString parFileOut = outFile + ".par.root";
+ TString geoFileOut = outFile + ".geo.root";
+
+ outFile += ".rec.root";
if ("" != sInpFile) inFile = sInpFile;
// ---------------------------------------------
@@ -206,7 +210,7 @@ Bool_t mcbm_event_reco_L1(UInt_t uRunId = 2570,
// Define output file for FairMonitor histograms
TString monitorFile{outFile};
- monitorFile.ReplaceAll("reco", "reco.monitor");
+ monitorFile.ReplaceAll(".rec.", ".rec.monitor.");
FairMonitor::GetMonitor()->EnableMonitor(kTRUE, monitorFile);
// ------------------------------------------------------------------------
@@ -656,14 +660,12 @@ Bool_t mcbm_event_reco_L1(UInt_t uRunId = 2570,
parIo3->open(parFileOut.Data(), "RECREATE");
// ------------------------------------------------------------------------
rtdb->setOutput(parIo3);
- rtdb->saveOutput();
- rtdb->print();
// ----- Run initialisation -------------------------------------------
std::cout << std::endl;
std::cout << "-I- " << myName << ": Initialise run" << std::endl;
run->Init();
-
+ rtdb->print();
// ------------------------------------------------------------------------
@@ -675,6 +677,11 @@ Bool_t mcbm_event_reco_L1(UInt_t uRunId = 2570,
// ----- Finish -------------------------------------------------------
+
+ rtdb->print();
+ rtdb->saveOutput();
+ run->CreateGeometryFile(geoFileOut);
+
timer.Stop();
FairMonitor::GetMonitor()->Print();
Double_t rtime = timer.RealTime();
diff --git a/macro/mcbm/mcbm_transport_boxgen.C b/macro/mcbm/mcbm_transport_boxgen.C
index b5b5ca88be025121ea45577b7d5356da4a413315..f47711dea95d1a2a8a48575d04db01c9e5d148c7 100644
--- a/macro/mcbm/mcbm_transport_boxgen.C
+++ b/macro/mcbm/mcbm_transport_boxgen.C
@@ -146,6 +146,12 @@ void mcbm_transport_boxgen(Int_t nEvents = 10,
run.SetParFileName(parFile);
run.SetGeoFileName(geoFile);
run.LoadSetup(setupName);
+
+ CbmSetup* geoSetup = CbmSetup::Instance();
+
+ // You can modify the pre-defined setup
+ // geoSetup->SetActive(ECbmModuleId::kMuch, kFALSE);
+
run.SetField(new CbmFieldConst());
run.SetTarget(targetElement, targetThickness, targetDiameter, targetPosX, targetPosY, targetPosZ,
targetRotY * TMath::DegToRad());
diff --git a/reco/alignment/CbmBbaAlignmentTask.cxx b/reco/alignment/CbmBbaAlignmentTask.cxx
index c3139a2c996c54a1da4ed88322e5ddaa2c57d483..43720e7f70bea7f646865e70d07d8fad15108b34 100644
--- a/reco/alignment/CbmBbaAlignmentTask.cxx
+++ b/reco/alignment/CbmBbaAlignmentTask.cxx
@@ -168,15 +168,15 @@ InitStatus CbmBbaAlignmentTask::Init()
ca::Framework* l1 = CbmL1::Instance()->fpAlgo;
const ca::Parameters<ca::fvec>& l1Param = l1->GetParameters();
- fNalignmentBodies = l1Param.GetNstationsActive();
+ fNtrackingStations = l1Param.GetNstationsActive();
TDirectory* curDirectory = gDirectory;
fHistoDir->cd();
- for (int i = -1; i < fNalignmentBodies; i++) {
- const char* n1 = Form("Body%d", i);
- const char* n2 = Form(", body %d", i);
+ for (int i = -1; i < fNtrackingStations; i++) {
+ const char* n1 = Form("Station%d", i);
+ const char* n2 = Form(", station %d", i);
if (i == -1) {
n1 = "All";
n2 = "";
@@ -283,7 +283,7 @@ void CbmBbaAlignmentTask::Exec(Option_t* /*opt*/)
break;
}
t.MakeConsistent();
- t.fUnalignedTrack.fMsQp0 = 1. / 0.1;
+ t.fUnalignedTrack.fMsQp0 = 1. / 0.5;
t.fUnalignedTrack.fIsMsQp0Set = true;
t.fUnalignedTrack.fNodes[t.fUnalignedTrack.fFirstHitNode].fTrack.Qp() = 0.;
@@ -324,13 +324,16 @@ void CbmBbaAlignmentTask::ApplyAlignment(const std::vector<double>& par)
for (unsigned int in = 0; in < t.fUnalignedTrack.fNodes.size(); in++) {
CbmKfTrackFitter::FitNode& node = t.fUnalignedTrack.fNodes[in];
CbmKfTrackFitter::FitNode& nodeAligned = t.fAlignedTrack.fNodes[in];
- int iAlignmentModule = node.fReference1;
- if (iAlignmentModule < 0) continue;
- assert(nodeAligned.fReference1 == iAlignmentModule);
+ int iSensor = node.fReference1;
+ assert(iSensor >= 0 && iSensor < (int) fSensors.size());
+ assert(nodeAligned.fReference1 == node.fReference1);
- double dx = par[3 * iAlignmentModule + 0];
- double dy = par[3 * iAlignmentModule + 1];
- double dz = par[3 * iAlignmentModule + 2];
+ Sensor& s = fSensors[iSensor];
+ if (s.fAlignmentBody < 0) continue;
+
+ double dx = par[3 * s.fAlignmentBody + 0];
+ double dy = par[3 * s.fAlignmentBody + 1];
+ double dz = par[3 * s.fAlignmentBody + 2];
nodeAligned.fMxy.SetX(node.fMxy.X() + ca::fvec(dx));
nodeAligned.fMxy.SetY(node.fMxy.Y() + ca::fvec(dy));
@@ -394,8 +397,8 @@ double CbmBbaAlignmentTask::CostFunction(const std::vector<double>& par)
if (fFixedNdf > 0 && fNdfTotal != fFixedNdf) {
cost = 1.e30;
}
- LOG(info) << "BBA: cost function: ndf " << fNdfTotal << ", cost " << cost
- << ", diff to ideal cost: " << cost - fCostIdeal;
+ LOG(info) << "BBA: calculate cost function: ndf " << fNdfTotal << ", cost " << cost
+ << ", diff to ideal cost: " << cost - fCostIdeal << ", diff to initial cost: " << cost - fCostInitial;
return cost;
}
@@ -408,15 +411,60 @@ void CbmBbaAlignmentTask::Finish()
LOG(info) << "BBA: start the alignment procedure with " << fTracks.size() << " tracks ...";
+ // collect the sensors from the hits
+ // TODO: read it from the detector setup
+
+
+ std::set<Sensor> sensorSet;
for (auto& t : fTracks) {
for (auto& n : t.fUnalignedTrack.fNodes) {
- // TODO: get the body index from n.fHitSystemId, n.fHitAddress
- int iBody = n.fMaterialLayer; // material layer is currently equal to the active tracking station
- n.fReference1 = iBody; // reference1 will be the body index
+ Sensor s;
+ s.fSystemId = n.fHitSystemId;
+ s.fSensorId = n.fHitAddress;
+ // TODO: get the station index from n.fHitSystemId, n.fHitAddress
+ s.fTrackingStation = n.fMaterialLayer;
+ sensorSet.insert(s);
+ }
+ }
+
+ fSensors.clear();
+ for (auto& s : sensorSet) { // convert set to a vector
+ fSensors.push_back(s);
+ }
+
+ for (auto& t : fTracks) {
+ for (auto& n : t.fUnalignedTrack.fNodes) {
+ Sensor s;
+ s.fSystemId = n.fHitSystemId;
+ s.fSensorId = n.fHitAddress;
+ auto iter = sensorSet.find(s);
+ assert(iter != sensorSet.end());
+ int iSensor = std::distance(sensorSet.begin(), iter);
+ n.fReference1 = iSensor;
}
t.fAlignedTrack = t.fUnalignedTrack;
}
+ if (1) { // one alignment body per tracking station
+
+ fNalignmentBodies = fNtrackingStations;
+
+ for (auto& s : fSensors) {
+ s.fAlignmentBody = s.fTrackingStation;
+ }
+ }
+ else { // one alignment body per sensor
+
+ fNalignmentBodies = fSensors.size();
+
+ for (int unsigned i = 0; i < fSensors.size(); i++) {
+ fSensors[i].fAlignmentBody = i;
+ }
+ }
+
+ LOG(info) << "BBA: " << fSensors.size() << " sensors, " << fNalignmentBodies << " alignment bodies";
+
+
int nParameters = 3 * fNalignmentBodies; // x, y, z
std::vector<bba::Parameter> par(nParameters);
@@ -425,7 +473,7 @@ void CbmBbaAlignmentTask::Finish()
bba::Parameter& p = par[ip];
p.SetActive(0);
p.SetValue(0.);
- p.SetBoundaryMin(-20.); // +-1 cm alignment range
+ p.SetBoundaryMin(-20.); // +-20 cm alignment range
p.SetBoundaryMax(20.);
p.SetStepMin(1.e-6);
p.SetStepMax(0.1);
@@ -434,16 +482,25 @@ void CbmBbaAlignmentTask::Finish()
// set active parameters
- for (int i = 1; i < fNalignmentBodies - 1; i++) {
- if (fTrackingMode == kSts && i == fNalignmentBodies / 2) {
- par[3 * i + 0].SetActive(0);
- par[3 * i + 1].SetActive(1);
- par[3 * i + 2].SetActive(0);
+ for (auto& s : fSensors) {
+ if (s.fAlignmentBody < 0) {
+ continue;
+ }
+ if (s.fTrackingStation == 0) { // the first station
+ continue;
+ }
+ if (s.fTrackingStation == fNtrackingStations - 1) { // the last station
+ continue;
+ }
+ if (fTrackingMode == kSts && s.fTrackingStation == fNtrackingStations / 2) { // station in the middle for STS mode
+ par[3 * s.fAlignmentBody + 0].SetActive(0);
+ par[3 * s.fAlignmentBody + 1].SetActive(1);
+ par[3 * s.fAlignmentBody + 2].SetActive(0);
}
else {
- par[3 * i + 0].SetActive(1);
- par[3 * i + 1].SetActive(1);
- par[3 * i + 2].SetActive(1);
+ par[3 * s.fAlignmentBody + 0].SetActive(1);
+ par[3 * s.fAlignmentBody + 1].SetActive(1);
+ par[3 * s.fAlignmentBody + 2].SetActive(1);
}
}
@@ -455,7 +512,7 @@ void CbmBbaAlignmentTask::Finish()
bba::Parameter& pz = par[3 * is + 2];
// +- d cm misalignment
- double d = 0.;
+ double d = fSimulatedMisalignmentRange;
if (px.IsActive()) {
px.SetValue(gRandom->Uniform(2. * d) - d);
}
@@ -511,20 +568,22 @@ void CbmBbaAlignmentTask::Finish()
if (!node.fIsFitted || !node.fIsXySet) continue;
+ Sensor& s = fSensors[node.fReference1];
+
double x_residual = node.fMxy.X()[0] - node.fTrack.X()[0]; // this should be the difference vector
double y_residual = node.fMxy.Y()[0] - node.fTrack.Y()[0];
hResidualsBeforeAlignmentX[0]->Fill(x_residual);
hResidualsBeforeAlignmentY[0]->Fill(y_residual);
- hResidualsBeforeAlignmentX[node.fReference1 + 1]->Fill(x_residual);
- hResidualsBeforeAlignmentY[node.fReference1 + 1]->Fill(y_residual);
+ hResidualsBeforeAlignmentX[s.fTrackingStation + 1]->Fill(x_residual);
+ hResidualsBeforeAlignmentY[s.fTrackingStation + 1]->Fill(y_residual);
double x_pull = x_residual / sqrt(node.fMxy.Dx2()[0] + node.fTrack.GetCovariance(0, 0)[0]);
double y_pull = y_residual / sqrt(node.fMxy.Dy2()[0] + node.fTrack.GetCovariance(1, 1)[0]);
hPullsBeforeAlignmentX[0]->Fill(x_pull);
hPullsBeforeAlignmentY[0]->Fill(y_pull);
- hPullsBeforeAlignmentX[node.fReference1 + 1]->Fill(x_pull);
- hPullsBeforeAlignmentY[node.fReference1 + 1]->Fill(y_pull);
+ hPullsBeforeAlignmentX[s.fTrackingStation + 1]->Fill(x_pull);
+ hPullsBeforeAlignmentY[s.fTrackingStation + 1]->Fill(y_pull);
}
}
}
@@ -550,14 +609,14 @@ void CbmBbaAlignmentTask::Finish()
LOG(info) << "\n Misaligned parameters: ";
for (int is = 0; is < fNalignmentBodies; is++) {
const std::vector<double>& r = parMisaligned;
- LOG(info) << "Module " << is << ": x " << r[3 * is + 0] << " y " << r[3 * is + 1] << " z " << r[3 * is + 2];
+ LOG(info) << "Body " << is << ": x " << r[3 * is + 0] << " y " << r[3 * is + 1] << " z " << r[3 * is + 2];
}
LOG(info) << "\n Alignment results: ";
for (int is = 0; is < fNalignmentBodies; is++) {
const std::vector<double>& r = alignment.getResult();
- LOG(info) << "Module " << is << ": x " << r[3 * is + 0] << " y " << r[3 * is + 1] << " z " << r[3 * is + 2];
+ LOG(info) << "Body " << is << ": x " << r[3 * is + 0] << " y " << r[3 * is + 1] << " z " << r[3 * is + 2];
}
LOG(info) << "\n";
@@ -572,20 +631,22 @@ void CbmBbaAlignmentTask::Finish()
if (!node.fIsFitted || !node.fIsXySet) continue;
+ Sensor& s = fSensors[node.fReference1];
+
double x_residual = node.fMxy.X()[0] - node.fTrack.X()[0];
double y_residual = node.fMxy.Y()[0] - node.fTrack.Y()[0];
hResidualsAfterAlignmentX[0]->Fill(x_residual);
hResidualsAfterAlignmentY[0]->Fill(y_residual);
- hResidualsAfterAlignmentX[node.fReference1 + 1]->Fill(x_residual);
- hResidualsAfterAlignmentY[node.fReference1 + 1]->Fill(y_residual);
+ hResidualsAfterAlignmentX[s.fTrackingStation + 1]->Fill(x_residual);
+ hResidualsAfterAlignmentY[s.fTrackingStation + 1]->Fill(y_residual);
double x_pull = x_residual / sqrt(node.fMxy.Dx2()[0] + node.fTrack.GetCovariance(0, 0)[0]);
double y_pull = y_residual / sqrt(node.fMxy.Dy2()[0] + node.fTrack.GetCovariance(1, 1)[0]);
hPullsAfterAlignmentX[0]->Fill(x_pull);
hPullsAfterAlignmentY[0]->Fill(y_pull);
- hPullsAfterAlignmentX[node.fReference1 + 1]->Fill(x_pull);
- hPullsAfterAlignmentY[node.fReference1 + 1]->Fill(y_pull);
+ hPullsAfterAlignmentX[s.fTrackingStation + 1]->Fill(x_pull);
+ hPullsAfterAlignmentY[s.fTrackingStation + 1]->Fill(y_pull);
}
}
diff --git a/reco/alignment/CbmBbaAlignmentTask.h b/reco/alignment/CbmBbaAlignmentTask.h
index c11571319a103e01a2c7667dbe04a764173c080e..af78145a7b0a5e0f4dfabea15d9557bba86c5ea5 100644
--- a/reco/alignment/CbmBbaAlignmentTask.h
+++ b/reco/alignment/CbmBbaAlignmentTask.h
@@ -13,6 +13,7 @@
#define CbmBbaAlignmentTask_H
+#include "CbmDefs.h"
#include "CbmKfTrackFitter.h"
#include "FairTask.h"
#include "TString.h"
@@ -51,6 +52,8 @@ class CbmBbaAlignmentTask : public FairTask {
void SetMcbmTrackingMode() { fTrackingMode = TrackingMode::kMcbm; }
void SetStsTrackingMode() { fTrackingMode = TrackingMode::kSts; }
+ void SetSimulatedMisalignmentRange(double range) { fSimulatedMisalignmentRange = range; }
+
public:
enum TrackingMode
{
@@ -75,6 +78,23 @@ class CbmBbaAlignmentTask : public FairTask {
void MakeConsistent();
};
+ struct Sensor {
+ ECbmModuleId fSystemId{0};
+ int fSensorId{-1};
+ int fTrackingStation{-1};
+ int fAlignmentBody{-1};
+ bool operator<(const Sensor& other) const
+ {
+ if (fSystemId < other.fSystemId) return true;
+ if (fSystemId > other.fSystemId) return false;
+ return (fSensorId < other.fSensorId);
+ };
+ bool operator==(const Sensor& other) const
+ {
+ return (fSystemId == other.fSystemId) && (fSensorId == other.fSensorId);
+ };
+ };
+
private:
const CbmBbaAlignmentTask& operator=(const CbmBbaAlignmentTask&);
CbmBbaAlignmentTask(const CbmBbaAlignmentTask&);
@@ -112,15 +132,20 @@ class CbmBbaAlignmentTask : public FairTask {
Int_t fMaxNtracks{100000};
+ int fNtrackingStations{0};
int fNalignmentBodies{0};
double fCostIdeal{1.e10};
double fCostInitial{0.};
+ double fSimulatedMisalignmentRange{0.}; // misalignment range for simulated misalignment
+
double fChi2Total{0.};
long fNdfTotal{0};
long fFixedNdf{-1};
+ std::vector<Sensor> fSensors;
+
//histograms
std::vector<TH1F*> hResidualsBeforeAlignmentX{};