BmnRoot
Loading...
Searching...
No Matches
run8_reco_bmn.cxx
Go to the documentation of this file.
1#if !defined(__CLING__) || defined(__MAKECLING__)
2// ROOT includes
3#include "TKey.h"
4#include "TList.h"
5#include "TPRegexp.h"
6#include "TROOT.h"
7#include "TRandom3.h"
8#include "TStopwatch.h"
9#include "TString.h"
10#include "TSystem.h"
11
12// Fair includes
13#include "FairCave.h"
14#include "FairField.h"
15#include "FairFileSource.h"
16#include "FairMagnet.h"
17#include "FairParAsciiFileIo.h"
18#include "FairParRootFileIo.h"
19#include "FairPipe.h"
20#include "FairRootFileSink.h"
21#include "FairRunAna.h"
22#include "FairRuntimeDb.h"
23#include "FairTask.h"
24#include "FairTrackParP.h"
25
26// BM@N includes
27#include "BmnBeamTracking.h"
28#include "BmnCSCHitMaker.h"
29#include "BmnCounterTask.h"
30#include "BmnDchHitProducer.h"
31#include "BmnDchTrackFinder.h"
32#include "BmnFieldMap.h"
33#include "BmnFileSource.h"
34#include "BmnFillDstTask.h"
35#include "BmnFunctionSet.h"
36#include "BmnGemResiduals.h"
37#include "BmnGemStripHitMaker.h"
38#include "BmnGemTrack.h"
39#include "BmnGlobalTracking.h"
40#include "BmnInnerTrackingRun7.h"
41#include "BmnLANDHitProducer.h"
42#include "BmnMwpcHitFinder.h"
43#include "BmnMwpcHitProducer.h"
44#include "BmnMwpcTrackFinder.h"
45#include "BmnNewFieldMap.h"
46#include "BmnPid.h"
47#include "BmnSiBTHitMaker.h"
48#include "BmnSiResiduals.h"
49#include "BmnSiliconHitMaker.h"
50#include "BmnTof1HitProducer.h"
51#include "BmnTofHitProducer.h"
52#include "BmnTriggersCheck.h"
53#include "BmnVertexFinder.h"
54#include "BmnZdcAnalyzer.h"
55#include "DstEventHeader.h"
56#include "MpdGetNumEvents.h"
57#include "UniRun.h"
58// L1 tracking
59#include "BmnToCbmHitConverter.h"
60#include "CbmFindPrimaryVertex.h"
61#include "CbmKF.h"
62#include "CbmL1.h"
63#include "CbmL1StsTrackFinder.h"
64#include "CbmPVFinderKF.h"
65#include "CbmStsFindTracks.h"
66#include "CbmStsTrackFinder.h"
67
68#include <iostream>
69using namespace std;
70#endif
71
72// -----------------------------------------------------------------------------
73// Macro for reconstruction of simulated or experimental events.
74//
75// inputFileName - input file with data.
76//
77// To process experimental data, you must use 'runN-NNN:'-like prefix
78// and then the geometry will be obtained from the Unified Condition Database.
79//
80// bmndstFileName - output file with reconstructed data.
81//
82// nStartEvent - number of first event to process (starts with zero), default: 0.
83//
84// nEvents - number of events to process, 0 - all events of given file will be
85// processed, default: 10000.
86//
87// alignCorrFileName - argument for choosing input file with the alignment
88// corrections.
89//
90// If alignCorrFileName == 'default', (case insensitive) then corrections are
91// retrieved from UniDb according to the running period and run number.
92//
93// If alignCorrFileName == '', then no corrections are applied at all.
94//
95// If alignCorrFileName == '<path>/<file-name>', then the corrections are taken
96// from that file.
97
98// #include "../run/bmnloadlibs.C"
99#define L1
100
101void run8_reco_bmn(TString inputFileName = "$VMCWORKDIR/macro/run8/bmnsim.root",
102 TString bmndstFileName = "$VMCWORKDIR/macro/run8/bmndst.root",
103 Int_t nStartEvent = 0,
104 Int_t nEvents = 10)
105{
106 gDebug = 0; // Debug option
107 // Verbosity level (0 = quiet (progress bar), 1 = event-level, 2 = track-level, 3 = full debug)
108 Int_t iVerbose = 0;
109
110 // ----- Timer ---------------------------------------------------------
111 TStopwatch timer;
112 timer.Start();
113
114 // check input file exists
115 if (!BmnFunctionSet::CheckFileExist(inputFileName, 1))
116 exit(-1);
117
118 // ----- Reconstruction run --------------------------------------------
119 FairRunAna* fRunAna = new FairRunAna();
120 fRunAna->SetEventHeader(new DstEventHeader());
121
122 Bool_t isField =
123 (inputFileName.Contains("noField")) ? kFALSE : kTRUE; // flag for tracking (to use mag.field or not)
124 Bool_t isTarget = kTRUE; // kTRUE; // flag for tracking (run with target or not)
125 Bool_t isExp = !BmnFunctionSet::isSimulationFile(
126 inputFileName); // flag for hit finder (to create digits or take them from data-file)
127
128 // Declare input source as simulation file or experimental data
129 FairSource* fFileSource;
130
131 // -1 means use of the BM@N-setup when processing MC-input
132 // DO NOT change it manually!
133 Int_t run_period = 8, run_number = -1;
134 Double_t fieldScale = 0.;
135 if (!isExp) // for simulation files
136 fFileSource = new FairFileSource(inputFileName);
137 else // for experimental files
138 {
139 // set source as raw root data file (without additional directories)
140 fFileSource = new BmnFileSource(inputFileName, run_period, run_number);
141
142 // get geometry for run
143 gRandom->SetSeed(0);
144 TString geoFileName = Form("current_geo_file_%d.root", UInt_t(gRandom->Integer(UINT32_MAX)));
145 Int_t res_code = UniRun::ReadGeometryFile(run_period, run_number, (char*)geoFileName.Data());
146 if (res_code != 0) {
147 cout << "ERROR: could not read geometry file from the database" << endl;
148 exit(-3);
149 }
150
151 // get gGeoManager from ROOT file (if required)
152 TFile* geoFile = new TFile(geoFileName, "READ");
153 if (!geoFile->IsOpen()) {
154 cout << "ERROR: could not open ROOT file with geometry: " + geoFileName << endl;
155 exit(-4);
156 }
157 TList* keyList = geoFile->GetListOfKeys();
158 TIter next(keyList);
159 TKey* key = (TKey*)next();
160 TString className(key->GetClassName());
161 if (className.BeginsWith("TGeoManager"))
162 key->ReadObj();
163 else {
164 cout << "ERROR: TGeoManager is not top element in geometry file " + geoFileName << endl;
165 exit(-5);
166 }
167
168 // set magnet field with factor corresponding to the given run
169 UniRun* pCurrentRun = UniRun::GetRun(run_period, run_number);
170 if (pCurrentRun == 0)
171 exit(-6);
172 Double_t* field_voltage = pCurrentRun->GetFieldVoltage();
173 if (field_voltage == NULL) {
174 cout << "ERROR: no field voltage was found for run " << run_period << ":" << run_number << endl;
175 exit(-7);
176 }
177 Double_t map_current = 55.87;
178 if (*field_voltage < 10) {
179 fieldScale = 0;
180 isField = kFALSE;
181 } else
182 fieldScale = (*field_voltage) / map_current;
183
184 BmnFieldMap* magField = new BmnNewFieldMap("field_sp41v5_ascii_Extrap.root");
185 magField->SetScale(fieldScale);
186 magField->Init();
187 fRunAna->SetField(magField);
188 isExp = kTRUE;
189 TString targ = "-";
190 if (pCurrentRun->GetTargetParticle() == "") {
191 isTarget = kFALSE;
192 } else {
193 targ = pCurrentRun->GetTargetParticle();
194 isTarget = kTRUE;
195 }
196 TString beam = pCurrentRun->GetBeamParticle();
197 if (beam == "")
198 beam = "-";
199 cout << "\n\n|||||||||||||||| EXPERIMENTAL RUN SUMMARY ||||||||||||||||" << endl;
200 cout << "||\t\t\t\t\t\t\t||" << endl;
201 cout << "||\t\tPeriod:\t\t" << run_period << "\t\t\t||" << endl;
202 cout << "||\t\tNumber:\t\t" << run_number << "\t\t\t||" << endl;
203 cout << "||\t\tBeam:\t\t" << beam << "\t\t\t||" << endl;
204 cout << "||\t\tTarget:\t\t" << targ << "\t\t\t||" << endl;
205 cout << "||\t\tField scale:\t" << setprecision(4) << fieldScale << "\t\t\t||" << endl;
206 cout << "||\t\t\t\t\t\t\t||" << endl;
207 cout << "||||||||||||||||||||||||||||||||||||||||||||||||||||||||||\n\n" << endl;
208 remove(geoFileName.Data());
209 }
210
211 fRunAna->SetSource(fFileSource);
212 // if directory for the output file does not exist, then create
213 if (BmnFunctionSet::CreateDirectoryTree(bmndstFileName, 1) < 0)
214 exit(-2);
215 fRunAna->SetSink(new FairRootFileSink(bmndstFileName));
216 fRunAna->SetGenerateRunInfo(false);
217
218 // if nEvents is equal 0 then all events of the given file starting with "nStartEvent" should be processed
219 if (nEvents == 0)
220 nEvents = MpdGetNumEvents::GetNumROOTEvents((char*)inputFileName.Data()) - nStartEvent;
221
222 // Digitisation files.
223 // Add TObjectString file names to a TList which is passed as input to the FairParAsciiFileIo.
224 // The FairParAsciiFileIo will create on fly a concatenated input parameter file, which is then used during the
225 // reconstruction.
226 TList* parFileNameList = new TList();
227
228 // ====================================================================== //
229 // === MWPC hit finder === //
230 // ====================================================================== //
231 // if(!isExp) {
232 // BmnMwpcHitProducer *mwpcHP = new BmnMwpcHitProducer();
233 // fRunAna->AddTask(mwpcHP);
234 // }
235 // BmnMwpcHitFinder* mwpcHM = new BmnMwpcHitFinder(isExp, run_period, run_number);
236 // fRunAna->AddTask(mwpcHM);
237
238 // ====================================================================== //
239 // === SiBT hit finder === //
240 // ====================================================================== //
241 BmnSiBTHitMaker* sibtHM = new BmnSiBTHitMaker(run_period, run_number, isExp);
242 fRunAna->AddTask(sibtHM);
243
244 // ====================================================================== //
245 // === Silicon hit finder === //
246 // ====================================================================== //
247 BmnSiliconHitMaker* siliconHM = new BmnSiliconHitMaker(run_period, run_number, isExp);
248 fRunAna->AddTask(siliconHM);
249
250 // ====================================================================== //
251 // === GEM hit finder === //
252 // ====================================================================== //
253 BmnGemStripHitMaker* gemHM = new BmnGemStripHitMaker(run_period, run_number, isExp);
254 gemHM->SetHitMatching(kTRUE);
255 fRunAna->AddTask(gemHM);
256
257 // ====================================================================== //
258 // === CSC hit finder === //
259 // ====================================================================== //
260 BmnCSCHitMaker* cscHM = new BmnCSCHitMaker(run_period, run_number, isExp);
261 if (!isExp)
262 cscHM->SetCurrentConfig(BmnCSCConfiguration::Run8); // set explicitly
263 cscHM->SetHitMatching(kTRUE);
264 fRunAna->AddTask(cscHM);
265
266 // ====================================================================== //
267 // === TOF1 hit finder === //
268 // ====================================================================== //
269 BmnTof1HitProducer* tof1HP = new BmnTof1HitProducer("TOF1", !isExp, iVerbose, kFALSE);
270 tof1HP->SetPeriodRun(run_period, run_number);
271 fRunAna->AddTask(tof1HP);
272
273 // ====================================================================== //
274 // === TOF2 hit finder === //
275 // ====================================================================== //
276 BmnTofHitProducer* tof2HP = new BmnTofHitProducer("TOF", "TOF700_geometry_run7.txt", !isExp, iVerbose, kFALSE);
277 tof2HP->SetTimeResolution(0.115);
278 tof2HP->SetProtonTimeCorrectionFile("bmn_run9687_digi_calibration.root");
279 tof2HP->SetMCTimeFile("TOF700_MC_argon_qgsm_time_run7.txt");
280 tof2HP->SetMainStripSelection(0); // 0 - minimal time, 1 - maximal amplitude
281 tof2HP->SetSelectXYCalibration(2); // 0 - Petukhov, 1 - Panin SRC, 2 - Petukhov Argon (default)
282 tof2HP->SetTimeMin(-2.f); // minimal digit time
283 tof2HP->SetTimeMax(+39.f); // Maximal digit time
284 tof2HP->SetDiffTimeMaxSmall(1.2f); // Abs maximal difference for small chambers
285 tof2HP->SetDiffTimeMaxBig(3.5f); // Abs maximal difference for big chambers
286 fRunAna->AddTask(tof2HP);
287
288 // ====================================================================== //
289 // === DCH hit finder === //
290 // ====================================================================== //
291 // if(!isExp) {
292 // BmnDchHitProducer *dchHP = new BmnDchHitProducer();
293 // fRunAna->AddTask(dchHP);
294 // }
295
296 // ====================================================================== //
297 // === ZDC === //
298 // ====================================================================== //
299 BmnZdcAnalyzer* zdcAna = new BmnZdcAnalyzer();
300 fRunAna->AddTask(zdcAna);
301
302#ifdef L1
303 FairTask* hitConverter = new BmnToCbmHitConverter(iVerbose);
304 fRunAna->AddTask(hitConverter);
305
306 // ====================================================================== //
307 // === STS track finding === //
308 // ====================================================================== //
309 CbmKF* kalman = new CbmKF();
310 fRunAna->AddTask(kalman);
311
312 CbmL1* l1 = new CbmL1();
313 TString stsMatBudgetFile =
314 ""; // paramDir + "/sts/sts_matbudget_v12b_12344444.root"; // paramDir + "/sts_matbudget_var_fr.root";
315 l1->SetMaterialBudgetFileName(stsMatBudgetFile);
316 fRunAna->AddTask(l1);
317
318 CbmStsTrackFinder* stsTrackFinder = new CbmL1StsTrackFinder();
319 FairTask* stsFindTracks = new CbmStsFindTracks(iVerbose, stsTrackFinder);
320 fRunAna->AddTask(stsFindTracks);
321#else
322 // ====================================================================== //
323 // === Tracking (InnerTracker) === //
324 // ====================================================================== //
325 BmnInnerTrackingRun7* innerTF = new BmnInnerTrackingRun7(run_number, isField, isTarget);
326 innerTF->SetFiltration(isExp); // we use filtration for experimental data only now
327 fRunAna->AddTask(innerTF);
328#endif
329
330 // ====================================================================== //
331 // === Tracking (BEAM) === //
332 // ====================================================================== //
333 BmnBeamTracking* beamTF = new BmnBeamTracking(run_period);
334 fRunAna->AddTask(beamTF);
335
336 // ====================================================================== //
337 // === Tracking (MWPC) === //
338 // ====================================================================== //
339 BmnMwpcTrackFinder* mwpcTF = new BmnMwpcTrackFinder(isExp, run_period, run_number);
340 fRunAna->AddTask(mwpcTF);
341
342 // ====================================================================== //
343 // === Tracking (DCH) === //
344 // ====================================================================== //
345 BmnDchTrackFinder* dchTF = new BmnDchTrackFinder(run_period, run_number, isExp);
346 dchTF->SetTransferFunction("transfer_func.txt");
347 fRunAna->AddTask(dchTF);
348
349 // ====================================================================== //
350 // === Global Tracking === //
351 // ====================================================================== //
352 Bool_t doAlign = kTRUE;
353 if (!isExp)
354 doAlign = kFALSE;
355 BmnGlobalTracking* glTF = new BmnGlobalTracking(isField, isExp, kFALSE /*doAlign*/);
356 fRunAna->AddTask(glTF);
357
358 // ====================================================================== //
359 // === Primary vertex finding === //
360 // ====================================================================== //
361#ifdef L1
362 CbmPrimaryVertexFinder* pvFinder = new CbmPVFinderKF();
363 CbmFindPrimaryVertex* findVertex = new CbmFindPrimaryVertex(pvFinder);
364 fRunAna->AddTask(findVertex);
365#else
366 BmnVertexFinder* gemVF = new BmnVertexFinder(run_period, isField);
367 fRunAna->AddTask(gemVF);
368#endif
369
370 // ====================================================================== //
371 // === PID procedure === //
372 // ====================================================================== //
373 BmnPid* pid_analyser = new BmnPid();
374 fRunAna->AddTask(pid_analyser);
375
376 // ====================================================================== //
377 // === Residual analysis === //
378 // ====================================================================== //
379 // BmnResiduals* res = new BmnResiduals(run_period, run_number, isField);
380 // fRunAna->AddTask(res);
381
382 // Fill DST Event Header (if iVerbose = 0, then print progress bar)
383 BmnFillDstTask* dst_task = new BmnFillDstTask(nEvents);
384 dst_task->SetRunNumber(run_period, run_number);
385 fRunAna->AddTask(dst_task);
386
387 // ----- Parameter database --------------------------------------------
388 FairRuntimeDb* rtdb = fRunAna->GetRuntimeDb();
389 FairParRootFileIo* parIo1 = new FairParRootFileIo();
390 FairParAsciiFileIo* parIo2 = new FairParAsciiFileIo();
391 parIo1->open(inputFileName.Data());
392 parIo2->open(parFileNameList, "in");
393 rtdb->setFirstInput(parIo1);
394 rtdb->setSecondInput(parIo2);
395 rtdb->setOutput(parIo1);
396 rtdb->saveOutput();
397 // -------------------------------------------------------------------------
398 // ----- Initialize and run --------------------------------------------
399 fRunAna->GetMainTask()->SetVerbose(iVerbose);
400 fRunAna->Init();
401 cout << "Starting run" << endl;
402 fRunAna->Run(nStartEvent, nStartEvent + nEvents);
403 // -------------------------------------------------------------------------
404 // ----- Finish --------------------------------------------------------
405 timer.Stop();
406 Double_t rtime = timer.RealTime();
407 Double_t ctime = timer.CpuTime();
408 cout << endl << endl;
409 cout << "Macro finished successfully." << endl; // marker of successful execution for CDASH
410 cout << "Input file is " + inputFileName << endl;
411 cout << "Output file is " + bmndstFileName << endl;
412 cout << "Real time " << rtime << " s, CPU time " << ctime << " s" << endl;
413 cout << endl;
414 // ------------------------------------------------------------------------
415 // delete fFileSource;
416 // delete fRunAna;
417}
418
419int main(int argc, char** arg)
420{
421 run8_reco_bmn();
422}
BmnCSCHitMaker::SetCurrentConfig
void SetCurrentConfig(BmnCSCConfiguration::CSC_CONFIG config)
Definition BmnCSCHitMaker.h:37
BmnCSCHitMaker::SetHitMatching
void SetHitMatching(Bool_t opt=kTRUE)
Definition BmnCSCHitMaker.h:35
BmnDchTrackFinder::SetTransferFunction
void SetTransferFunction(TString name)
Definition BmnDchTrackFinder.h:54
BmnFieldMap::SetScale
void SetScale(Double_t factor)
Definition BmnFieldMap.h:55
BmnFieldMap::Init
virtual void Init()
Definition BmnFieldMap.cxx:93
BmnFillDstTask::SetRunNumber
void SetRunNumber(Int_t period_number, Int_t run_number)
Definition BmnFillDstTask.h:58
BmnFunctionSet::CreateDirectoryTree
static int CreateDirectoryTree(TString &fileName, int iVerbose=0, EAccessMode mode=kWritePermission)
Definition BmnFunctionSet.cxx:56
BmnFunctionSet::isSimulationFile
static int isSimulationFile(TString fileName="")
Definition BmnFunctionSet.cxx:88
BmnFunctionSet::CheckFileExist
static int CheckFileExist(TString &fileName, int iVerbose=0, EAccessMode mode=kFileExists)
Definition BmnFunctionSet.cxx:16
BmnGemStripHitMaker::SetHitMatching
void SetHitMatching(Bool_t opt=kTRUE)
Definition BmnGemStripHitMaker.h:36
BmnTof1HitProducer::SetPeriodRun
void SetPeriodRun(Int_t p, Int_t r)
Definition BmnTof1HitProducer.h:83
BmnTofHitProducer::SetTimeResolution
void SetTimeResolution(Double_t sigma)
Definition BmnTofHitProducer.h:88
BmnTofHitProducer::SetDiffTimeMaxSmall
void SetDiffTimeMaxSmall(Double_t ds=0.)
Definition BmnTofHitProducer.h:98
BmnTofHitProducer::SetTimeMax
void SetTimeMax(Double_t ma=+15.)
Definition BmnTofHitProducer.h:97
BmnTofHitProducer::SetSelectXYCalibration
void SetSelectXYCalibration(Int_t s=0)
Definition BmnTofHitProducer.h:95
BmnTofHitProducer::SetMCTimeFile
void SetMCTimeFile(const char *file)
Definition BmnTofHitProducer.h:113
BmnTofHitProducer::SetDiffTimeMaxBig
void SetDiffTimeMaxBig(Double_t db=0.)
Definition BmnTofHitProducer.h:99
BmnTofHitProducer::SetMainStripSelection
void SetMainStripSelection(Int_t s=0)
Definition BmnTofHitProducer.h:94
BmnTofHitProducer::SetProtonTimeCorrectionFile
void SetProtonTimeCorrectionFile(const char *file)
Definition BmnTofHitProducer.h:114
BmnTofHitProducer::SetTimeMin
void SetTimeMin(Double_t mi=-2.)
Definition BmnTofHitProducer.h:96
CbmKF
Definition CbmKF.h:28
CbmL1
Definition CbmL1.h:49
CbmL1::SetMaterialBudgetFileName
void SetMaterialBudgetFileName(TString s)
Definition CbmL1.h:78
MpdGetNumEvents::GetNumROOTEvents
static Int_t GetNumROOTEvents(const char *filename, const char *treename="", int iVerbose=2)
Definition MpdGetNumEvents.cxx:511
UniRun::GetRun
static UniRun * GetRun(int period_number, int run_number)
get run from the database
Definition UniRun.cxx:177
UniRun::ReadGeometryFile
static int ReadGeometryFile(int period_number, int run_number, char *geo_file_path, bool usePrevGeometryIfMissing=false)
Definition UniRun.cxx:1105
UniRun::GetBeamParticle
TString GetBeamParticle()
get beam particle of the current run
Definition UniRun.h:113
UniRun::GetFieldVoltage
double * GetFieldVoltage()
get field voltage of the current run
Definition UniRun.h:125
UniRun::GetTargetParticle
TString GetTargetParticle()
get target particle of the current run
Definition UniRun.h:115
std
STL namespace.
run8_reco_bmn
void run8_reco_bmn(TString inputFileName="$VMCWORKDIR/macro/run8/bmnsim.root", TString bmndstFileName="$VMCWORKDIR/macro/run8/bmndst.root", Int_t nStartEvent=0, Int_t nEvents=10)
Definition run8_reco_bmn.cxx:101
main
int main()
Definition run_reco_src.cxx:378