// -------------------------------------------------------------------------- // // Macro for reconstruction of simulated events with standard settings // // HitProducers in MVD, RICH, TRD, TOF, ECAL // Digitizer and HitFinder in STS // FAST MC for ECAL // STS track finding and fitting (L1 / KF) // TRD track finding and fitting (L1 / KF) // RICH ring finding (ideal) and fitting // Global track finding (ideal), rich assignment // Primary vertex finding (ideal) // Matching of reconstructed and MC tracks in STS, RICH and TRD // // V. Friese 24/02/2006 // Version 04/03/2015 (V. Friese) // // -------------------------------------------------------------------------- void mcbm_reco_nh(Int_t nEvents = 1000, TString cSys="lam", TString cEbeam="2.5gev", TString cCentr="-", Int_t iRun=0, const char* setup = "sis18_mcbm_20deg_long") { // ======================================================================== // Adjust this part according to your requirements // Verbosity level (0=quiet, 1=event level, 2=track level, 3=debug) Int_t iVerbose = 0; TString outDir = "data/"; TString inFile = outDir + setup + "_" + cSys + "." + cEbeam + "." + cCentr + ".mc." + Form("%05d",iRun) + ".root"; // Input file (MC events) TString parFile = outDir + setup + "_" + cSys + "." + cEbeam + "." + cCentr + ".params." + Form("%05d",iRun) + ".root"; // Parameter file TString outFile = outDir + setup + "_" + cSys + "." + cEbeam + "." + cCentr + ".eds." + Form("%05d",iRun) + ".root"; // Output file FairLogger::GetLogger()->SetLogScreenLevel("WARNING"); //FairLogger::GetLogger()->SetLogScreenLevel("INFO"); //FairLogger::GetLogger()->SetLogScreenLevel("DEBUG"); FairLogger::GetLogger()->SetLogVerbosityLevel("MEDIUM"); // Digitisation files. // Add TObjectString containing the different file names to // a TList which is passed as input to the FairParAsciiFileIo. // The FairParAsciiFileIo will take care to create on the fly // a concatenated input parameter file which is then used during // the reconstruction. TList *parFileList = new TList(); TString inDir = gSystem->Getenv("VMCWORKDIR"); TString paramDir = inDir + "/parameters/"; TString setupFile = inDir + "/geometry/setup/legacy/" + setup + "_setup.C"; TString setupFunct = setup; setupFunct += "_setup()"; gROOT->LoadMacro(setupFile); gInterpreter->ProcessLine(setupFunct); /* // ----- Load the geometry setup ------------------------------------- std::cout << std::endl; TString setupFile = inDir + "/geometry/setup/setup_" + setup+ ".C"; TString setupFunct = "setup_"; setupFunct = setupFunct + setup + "()"; std::cout << "-I- mcbm_reco: Loading macro " << setupFile << std::endl; gROOT->LoadMacro(setupFile); gROOT->ProcessLine(setupFunct); */ // --- STS digipar file is there only for L1. It is no longer required // --- for STS digitisation and should be eventually removed. //TObjString stsDigiFile = paramDir + stsDigi; //parFileList->Add(&stsDigiFile); //cout << "macro/run/run_reco.C using: " << stsDigi << endl; TObjString* trdDigiFile = new TObjString(paramDir + trdDigi); parFileList->Add(trdDigiFile); // TObjString trdDigiFile(paramDir + trdDigi); // parFileList->Add(&trdDigiFile); cout << "macro/run/run_reco.C using: " << trdDigi << endl; TObjString* tofDigiFile = new TObjString(paramDir + tofDigi); parFileList->Add(tofDigiFile); cout << "macro/mcbm/mcbm_reco.C using: " << tofDigi << endl; TObjString* tofDigiBdfFile = new TObjString(paramDir + tofDigiBdf); parFileList->Add(tofDigiBdfFile); cout << "macro/mcbm/mcbm_reco.C using: " << paramDir << tofDigiBdf << endl; // TObjString tofDigiFile = paramDir + tofDigi; // parFileList->Add(&tofDigiFile); // In general, the following parts need not be touched // ======================================================================== // ---- Debug option ------------------------------------------------- gDebug = 0; // ------------------------------------------------------------------------ // ----- Timer -------------------------------------------------------- TStopwatch timer; timer.Start(); // ------------------------------------------------------------------------ // ----- Reconstruction run ------------------------------------------- FairRunAna *run = new FairRunAna(); run->SetInputFile(inFile); run->SetOutputFile(outFile); // ------------------------------------------------------------------------ // ----- MC Data Manager ------------------------------------------------ CbmMCDataManager* mcManager=new CbmMCDataManager("MCManager", 1); mcManager->AddFile(inFile); run->AddTask(mcManager); // ------------------------------------------------------------------------ // ========================================================================= // === Detector Response Simulation (Digitiser) === // === (where available) === // ========================================================================= // ----- MVD Digitiser ------------------------------------------------- CbmMvdDigitizer* mvdDigitise = new CbmMvdDigitizer("MVD Digitiser", 0, iVerbose); //run->AddTask(mvdDigitise); // ------------------------------------------------------------------------- // ----- MVD Clusterfinder --------------------------------------------- CbmMvdClusterfinder* mvdCluster = new CbmMvdClusterfinder("MVD Clusterfinder", 0, iVerbose); //run->AddTask(mvdCluster); // ------------------------------------------------------------------------- // ----- STS digitizer ------------------------------------------------- // ----- The parameters of the STS digitizer are set such as to match // ----- those in the old digitizer. Change them only if you know what you // ----- are doing. Double_t dynRange = 40960.; // Dynamic range [e] Double_t threshold = 4000.; // Digitisation threshold [e] Int_t nAdc = 4096; // Number of ADC channels (12 bit) Double_t timeResolution = 5.; // time resolution [ns] Double_t deadTime = 9999999.; // infinite dead time (integrate entire event) Double_t noise = 0.; // ENC [e] Int_t digiModel = 1; // Model: 1 = uniform charge distribution along track CbmStsDigitize* stsDigi = new CbmStsDigitize(digiModel); stsDigi->SetParameters(dynRange, threshold, nAdc, timeResolution, deadTime, noise); run->AddTask(stsDigi); // ------------------------------------------------------------------------- // ========================================================================= // === MVD local reconstruction === // ========================================================================= // ----- MVD Hit Finder ------------------------------------------------ CbmMvdHitfinder* mvdHitfinder = new CbmMvdHitfinder("MVD Hit Finder", 0, iVerbose); mvdHitfinder->UseClusterfinder(kTRUE); //run->AddTask(mvdHitfinder); // ------------------------------------------------------------------------- // === End of MVD local reconstruction === // ========================================================================= // ========================================================================= // === STS local reconstruction === // ========================================================================= // ----- STS Cluster Finder -------------------------------------------- FairTask* stsClusterFinder = new CbmStsFindClusters(); run->AddTask(stsClusterFinder); // ------------------------------------------------------------------------- // ----- STS hit finder ------------------------------------------------ FairTask* stsFindHits = new CbmStsFindHits(); run->AddTask(stsFindHits); // ------------------------------------------------------------------------- // ----- STS hit matching ----------------------------------------------- // FairTask* stsMatchHits = new CbmStsMatchHits(); // run->AddTask(stsMatchHits); // ------------------------------------------------------------------------- // --- STS track finding ------------------------------------------------ CbmKF* kalman = new CbmKF(); run->AddTask(kalman); CbmL1* l1 = new CbmL1(); TString mvdMatBudgetFileName = paramDir + mvdMatBudget; TString stsMatBudgetFileName = paramDir + stsMatBudget; // l1->SetStsMaterialBudgetFileName(stsMatBudgetFileName.Data()); // l1->SetMvdMaterialBudgetFileName(mvdMatBudgetFileName.Data()); // run->AddTask(l1); CbmStsTrackFinder* stsTrackFinder = new CbmL1StsTrackFinder(); FairTask* stsFindTracks = new CbmStsFindTracks(iVerbose, stsTrackFinder); // run->AddTask(stsFindTracks); // ------------------------------------------------------------------------- // --- STS track matching ---------------------------------------------- // FairTask* stsMatchTracks = new CbmStsMatchTracks(iVerbose); // run->AddTask(stsMatchTracks); // ------------------------------------------------------------------------- // --- STS track fitting ----------------------------------------------- // CbmStsTrackFitter* stsTrackFitter = new CbmStsKFTrackFitter(); // FairTask* stsFitTracks = new CbmStsFitTracks(stsTrackFitter, iVerbose); // run->AddTask(stsFitTracks); // ------------------------------------------------------------------------- // === End of STS local reconstruction === // ========================================================================= // ========================================================================= // === TRD local reconstruction === // ========================================================================= CbmTrdRadiator *radiator = new CbmTrdRadiator(kTRUE,"K++"); FairTask* trdDigi = new CbmTrdDigitizerPRF(radiator); //run->AddTask(trdDigi); Double_t triggerThreshold = 0.5e-6; // SIS100 Bool_t triangularPads = false; // Bucharest triangular pad-plane layout CbmTrdClusterFinderFast* trdCluster = new CbmTrdClusterFinderFast(); trdCluster->SetNeighbourTrigger(true); trdCluster->SetTriggerThreshold(triggerThreshold); trdCluster->SetNeighbourRowTrigger(false); trdCluster->SetPrimaryClusterRowMerger(true); trdCluster->SetTriangularPads(triangularPads); //run->AddTask(trdCluster); CbmTrdHitProducerCluster* trdHit = new CbmTrdHitProducerCluster(); trdHit->SetTriangularPads(triangularPads); //run->AddTask(trdHit); /* Bool_t simpleTR = kTRUE; // use fast and simple version for TR production CbmTrdRadiator *radiator = new CbmTrdRadiator(simpleTR,"K++"); //"K++" : micro structured POKALON //"H++" : PE foam foils //"G30" : ALICE fibers 30 layers Bool_t triangularPads = false;// Bucharest triangular pad-plane layout //Double_t triggerThreshold = 0.5e-6;//SIS100 Double_t triggerThreshold = 1.0e-6;//SIS300 Double_t trdNoiseSigma_keV = 0.1; //default best matching to test beam PRF CbmTrdDigitizerPRF* trdDigiPrf = new CbmTrdDigitizerPRF(radiator); trdDigiPrf->SetTriangularPads(triangularPads); trdDigiPrf->SetNoiseLevel(trdNoiseSigma_keV); run->AddTask(trdDigiPrf); CbmTrdClusterFinderFast* trdCluster = new CbmTrdClusterFinderFast(); trdCluster->SetNeighbourTrigger(true); trdCluster->SetTriggerThreshold(triggerThreshold); trdCluster->SetNeighbourRowTrigger(false); trdCluster->SetPrimaryClusterRowMerger(true); trdCluster->SetTriangularPads(triangularPads); run->AddTask(trdCluster); CbmTrdHitProducerCluster* trdHit = new CbmTrdHitProducerCluster(); trdHit->SetTriangularPads(triangularPads); // run->AddTask(trdHit); */ // ------------------------------------------------------------------------- // === End of TRD local reconstruction === // ========================================================================= // ========================================================================= // === TOF local reconstruction === // ========================================================================= CbmTofDigitizerBDF* tofDigi = new CbmTofDigitizerBDF("TOF Digitizer BDF",iVerbose); tofDigi->SetOutputBranchPersistent("TofDigi", kTRUE); tofDigi->SetOutputBranchPersistent("TofDigiMatchPoints", kTRUE); tofDigi->SetInputFileName( paramDir + "/tof/test_bdf_input.root"); // tofDigi->SetHistoFileName( digiOutFile ); // Uncomment to save control histograms run->AddTask(tofDigi); CbmTofSimpClusterizer* tofCluster = new CbmTofSimpClusterizer("TOF Simple Clusterizer", 0); tofCluster->SetOutputBranchPersistent("TofHit", kTRUE); tofCluster->SetOutputBranchPersistent("TofDigiMatch", kTRUE); run->AddTask(tofCluster); // ------ TOF hit producer --------------------------------------------- // CbmTofHitProducerNew* tofHitProd = new CbmTofHitProducerNew("TOF HitProducerNew",iVerbose); // tofHitProd->SetInitFromAscii(kFALSE); // run->AddTask(tofHitProd); // ------------------------------------------------------------------------- // === End of TOF local reconstruction === // ========================================================================= // ========================================================================= // === Global tracking === // ========================================================================= CbmLitFindGlobalTracks* finder = new CbmLitFindGlobalTracks(); // Tracking method to be used // "branch" - branching tracking // "nn" - nearest neighbor tracking // "weight" - weighting tracking finder->SetTrackingType("branch"); // Hit-to-track merger method to be used // "nearest_hit" - assigns nearest hit to the track finder->SetMergerType("nearest_hit"); // run->AddTask(finder); // ----- Primary vertex finding --------------------------------------- CbmPrimaryVertexFinder* pvFinder = new CbmPVFinderKF(); CbmFindPrimaryVertex* findVertex = new CbmFindPrimaryVertex(pvFinder); //run->AddTask(findVertex); // ------------------------------------------------------------------------ // === End of global tracking === // ========================================================================= // ----------- TRD track Pid Wkn ---------------------- CbmTrdSetTracksPidWkn* trdSetTracksPidTask = new CbmTrdSetTracksPidWkn( "trdFindTracks", "trdFindTracks"); // run->AddTask(trdSetTracksPidTask); // ---------------------------------------------------- // ----------- TRD track Pid Ann ---------------------- CbmTrdSetTracksPidANN* trdSetTracksPidAnnTask = new CbmTrdSetTracksPidANN( "Ann", "Ann"); // run->AddTask(trdSetTracksPidAnnTask); // ---------------------------------------------------- // ----------- TRD track Pid Like ---------------------- // Since in the newest version of this method depends on the global // track the task has to move after the global tracking // FU 08.02.12 Switch the task off since the input file needed for the new geometry has to be generated first. // CbmTrdSetTracksPidLike* trdSetTracksPidLikeTask = // new CbmTrdSetTracksPidLike("Likelihood", "Likelihood"); // run->AddTask(trdSetTracksPidLikeTask); // ---------------------------------------------------- // ========================================================================= // === RICH reconstruction === // ========================================================================= if (richGeom.Length() != 0) // if RICH is defined { // ---------------------RICH Hit Producer ---------------------------------- CbmRichHitProducer* richHitProd = new CbmRichHitProducer(); run->AddTask(richHitProd); /* CbmRichHitProducer* richHitProd = new CbmRichHitProducer(); richHitProd->SetDetectorType(4); richHitProd->SetNofNoiseHits(220); richHitProd->SetCollectionEfficiency(1.0); richHitProd->SetSigmaMirror(0.06); run->AddTask(richHitProd); */ //-------------------------------------------------------------------------- //--------------------- RICH Reconstruction ---------------------------------- CbmRichReconstruction* richReco = new CbmRichReconstruction(); run->AddTask(richReco); // ------------------- RICH Ring matching --------------------------------- // CbmRichMatchRings* matchRings = new CbmRichMatchRings(); // run->AddTask(matchRings); // ------------------------------------------------------------------------- } // === End of RICH local reconstruction === // ========================================================================= /* // ========================================================================= // === ECAL reconstruction === // ========================================================================= // ----- ECAL hit producer ---------------------------------------------- CbmEcalHitProducerFastMC* ecalHitProd = new CbmEcalHitProducerFastMC( "ECAL Hitproducer"); run->AddTask(ecalHitProd); // ------------------------------------------------------------------------- // === End of ECAL reconstruction === // ========================================================================= */ // ========================================================================= // === Matching to Monte-carlo === // ========================================================================= // CbmMatchRecoToMC* matchTask = new CbmMatchRecoToMC(); // run->AddTask(matchTask); // === End of matching to Monte-Carlo === // ========================================================================= // ----- Parameter database -------------------------------------------- FairRuntimeDb* rtdb = run->GetRuntimeDb(); FairParRootFileIo* parIo1 = new FairParRootFileIo(); FairParAsciiFileIo* parIo2 = new FairParAsciiFileIo(); parIo1->open(parFile.Data()); parIo2->open(parFileList, "in"); rtdb->setFirstInput(parIo1); rtdb->setSecondInput(parIo2); rtdb->setOutput(parIo1); rtdb->saveOutput(); // ------------------------------------------------------------------------ // ----- Intialise and run -------------------------------------------- run->Init(); cout << "Starting run" << endl; run->Run(0, nEvents); // ------------------------------------------------------------------------ // ----- Finish ------------------------------------------------------- timer.Stop(); Double_t rtime = timer.RealTime(); Double_t ctime = timer.CpuTime(); cout << endl << endl; cout << "Macro finished succesfully." << endl; cout << "Output file is " << outFile << endl; cout << "Parameter file is " << parFile << endl; cout << "Real time " << rtime << " s, CPU time " << ctime << " s" << endl; cout << endl; // ------------------------------------------------------------------------ // delete run; cout << " Test passed" << endl; cout << " All ok " << endl; RemoveGeoManager(); }