// -------------------------------------------------------------------------- // // 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 // // -------------------------------------------------------------------------- void run_reco(Int_t nEvents = 1) { // ======================================================================== // Adjust this part according to your requirements // Verbosity level (0=quiet, 1=event level, 2=track level, 3=debug) Int_t iVerbose = 0; // Input file (MC events) TString inFile = "data/test.mc.root"; // Parameter file TString parFile = "data/params.root"; // STS digitisation file TString stsDigiFile = "sts_standard.gsi.digi.par"; // Output file TString outFile = "data/test.reco.root"; // In general, the following parts need not be touched // ======================================================================== // ---- Debug option ------------------------------------------------- gDebug = 0; // ------------------------------------------------------------------------ // ----- Timer -------------------------------------------------------- TStopwatch timer; timer.Start(); // ------------------------------------------------------------------------ // ---- Load libraries ------------------------------------------------- gROOT->LoadMacro("$VMCWORKDIR/gconfig/basiclibs.C"); basiclibs(); gSystem->Load("libGeoBase"); gSystem->Load("libParBase"); gSystem->Load("libBase"); gSystem->Load("libMCStack"); gSystem->Load("libField"); gSystem->Load("libGen"); gSystem->Load("libPassive"); gSystem->Load("libMvd"); gSystem->Load("libSts"); gSystem->Load("libRich"); gSystem->Load("libTrd"); gSystem->Load("libTof"); gSystem->Load("libEcal"); gSystem->Load("libGlobal"); gSystem->Load("libKF"); gSystem->Load("libL1"); // ------------------------------------------------------------------------ // ----- Reconstruction run ------------------------------------------- CbmRunAna *fRun= new CbmRunAna(); fRun->SetInputFile(inFile); fRun->SetOutputFile(outFile); // ------------------------------------------------------------------------ // ========================================================================= // === MVD local reconstruction === // ========================================================================= // ----- MVD Hitproducer ----------------------------------------------- CbmMvdHitProducer* hitProd = new CbmMvdHitProducer("MVDHitProducer", 0, iVerbose); fRun->AddTask(hitProd); // ------------------------------------------------------------------------- // === End of MVD local reconstruction === // ========================================================================= // ========================================================================= // === STS local reconstruction === // ========================================================================= // ----- STS digitizer ------------------------------------------------- CbmTask* stsDigitize = new CbmStsDigitize(iVerbose); fRun->AddTask(stsDigitize); // ------------------------------------------------------------------------- // ----- STS hit finding ------------------------------------------------ CbmTask* findHits = new CbmStsFindHits(iVerbose); fRun->AddTask(findHits); // ------------------------------------------------------------------------- // ----- STS hit matching ----------------------------------------------- CbmStsMatchHits* matchHits = new CbmStsMatchHits(iVerbose); fRun->AddTask(matchHits); // ------------------------------------------------------------------------- // --- STS track finding ------------------------------------------------ CbmKF* kalman = new CbmKF(); fRun->AddTask(kalman); CbmL1* l1 = new CbmL1(); fRun->AddTask(l1); CbmStsTrackFinder* trackFinder = new CbmL1StsTrackFinder(); CbmTask* findTracks = new CbmStsFindTracks(trackFinder, iVerbose); fRun->AddTask(findTracks); // ------------------------------------------------------------------------- // --- STS track matching ---------------------------------------------- CbmTask* matchTracks = new CbmStsMatchTracks(iVerbose); fRun->AddTask(matchTracks); // ------------------------------------------------------------------------- // --- STS track fitting ----------------------------------------------- CbmStsTrackFitter* trackFitter = new CbmStsKFTrackFitter(); CbmTask* fitTracks = new CbmStsFitTracks(trackFitter, iVerbose); fRun->AddTask(fitTracks); // ------------------------------------------------------------------------- // === End of STS local reconstruction === // ========================================================================= // ========================================================================= // === TRD local reconstruction === // ========================================================================= // ----- TRD hit producer ---------------------------------------------- Double_t trdSigmaX[] = {300, 400, 500}; // Resolution in x [mum] // Resolutions in y - station and angle dependent [mum] Double_t trdSigmaY1[] = {2700, 3700, 15000, 27600, 33000, 33000, 33000 }; Double_t trdSigmaY2[] = {6300, 8300, 33000, 33000, 33000, 33000, 33000 }; Double_t trdSigmaY3[] = {10300, 15000, 33000, 33000, 33000, 33000, 33000 }; Int_t trdNFoils = 100; // number of polyetylene foils Float_t trdDFoils = 0.001; // thickness of 1 foil [cm] Float_t trdDGap = 0.008; // thickness of gap between foils [cm] CbmTrdHitProducer* trdHitProd = new CbmTrdHitProducer("TRD Hitproducer", "TRD task"); trdHitProd->SetPar(trdNFoils, trdDFoils, trdDGap); trdHitProd->SetSigmaX(trdSigmaX); trdHitProd->SetSigmaY(trdSigmaY1, trdSigmaY2, trdSigmaY3); fRun->AddTask(trdHitProd); // ------------------------------------------------------------------------- // ----- TRD track finding --------------------------------------------- CbmTrdTrackFinder* trdTrackFinder = new CbmL1TrdTrackFinderSts(); ((CbmL1TrdTrackFinderSts*)trdTrackFinder)->SetVerbose(iVerbose); CbmTrdFindTracks* trdFindTracks = new CbmTrdFindTracks("TRD Track Finder"); trdFindTracks->UseFinder(trdTrackFinder); fRun->AddTask(trdFindTracks); // ------------------------------------------------------------------------- // ----- TRD track fitting --------------------------------------------- CbmTrdTrackFitter* trdTrackFitter = new CbmTrdTrackFitterKF(); ((CbmTrdTrackFitterKF*)trdTrackFitter)->SetVerbose(iVerbose); ((CbmTrdTrackFitterKF*)trdTrackFitter)->SetPid(211); CbmTrdFitTracks* trdFitTracks = new CbmTrdFitTracks("TRD track fitter", "TRD", trdTrackFitter); fRun->AddTask(trdFitTracks); // ------------------------------------------------------------------------- // ----- TRD track matching -------------------------------------------- CbmTrdMatchTracks* trdMatchTracks = new CbmTrdMatchTracks(iVerbose); fRun->AddTask(trdMatchTracks); // ------------------------------------------------------------------------- // === End of TRD local reconstruction === // ========================================================================= // ========================================================================= // === TOF local reconstruction === // ========================================================================= // ------ TOF hit producer --------------------------------------------- CbmTofHitProducer* tofHitProd = new CbmTofHitProducer("TOF HitProducer", iVerbose); fRun->AddTask(tofHitProd); // ------------------------------------------------------------------------- // === End of TOF local reconstruction === // ========================================================================= // ========================================================================= // === Global tracking === // ========================================================================= // ----- STS-TRD-TOF track merging ------------------------------------- CbmL1TrackMerger* trackMerger = new CbmL1TrackMerger(); trackMerger->SetMethod(1); CbmL1TofMerger* tofMerger = new CbmL1TofMerger(); CbmFindGlobalTracks* findGlobal = new CbmFindGlobalTracks(trackMerger, NULL, tofMerger, iVerbose); fRun->AddTask(findGlobal); // ------------------------------------------------------------------------ // ----- Primary vertex finding --------------------------------------- CbmPrimaryVertexFinder* pvFinder = new CbmPVFinderKF(); CbmFindPrimaryVertex* findVertex = new CbmFindPrimaryVertex(pvFinder); fRun->AddTask(findVertex); // ------------------------------------------------------------------------ // === End of global tracking === // ========================================================================= // ========================================================================= // === RICH reconstruction === // ========================================================================= // ---------------------RICH Hit Producer ---------------------------------- Double_t richPmtRad = 0.4; // PMT radius [cm] Double_t richPmtDist = 0.; // Distance between PMTs [cm] Int_t richDetType = 1; // Detector type Protvino Int_t richNoise = 220; // Number of noise points per event CbmRichHitProducer* richHitProd = new CbmRichHitProducer(richPmtRad, richPmtDist, richDetType, richNoise, iVerbose); fRun->AddTask(richHitProd); //-------------------------------------------------------------------------- //----------------------RICH Track Extrapolation --------------------------- Int_t richNSts = 4; // minimum number of STS hits for extrapolation Double_t richZPos = 300.; // z position for extrapolation [cm] CbmRichTrackExtrapolation* richExtra = new CbmRichTrackExtrapolationKF(richNSts, iVerbose); CbmRichExtrapolateTracks* richExtrapolate = new CbmRichExtrapolateTracks(); richExtrapolate->UseExtrapolation(richExtra,richZPos); fRun->AddTask(richExtrapolate); //-------------------------------------------------------------------------- //--------------------- Rich Track Projection to photodetector ------------- Int_t richZFlag = 1; // Projetion from IM plane (default) CbmRichProjectionProducer* richProj = new CbmRichProjectionProducer(iVerbose, richZFlag); fRun->AddTask(richProj); //-------------------------------------------------------------------------- //--------------------- RICH Ring Finding ---------------------------------- Int_t NIterations = 3; Int_t NBinsX[] = {250,250,150}; Int_t NBinsY[] = {100,100,60}; Int_t HTNBinsX[] = {150,150,150}; Int_t HTNBinsY[] = {60,60,60}; Int_t HTCut[] = {300,50,30}; Int_t HitCut[] = {15,3,3}; Double_t MaxDistance[] = {13.,13.,13.0}; Double_t MinDistance[] = {1.,1.,1.}; Double_t MinRadius[] = {6.,6.,3.0}; Double_t MaxRadius[] = {6.5,6.5,6.2}; Int_t NParts[] = {3,2,2}; Int_t SqLen = 1; CbmRichRingFinderHough* richFinder = new CbmRichRingFinderHough(iVerbose); richFinder->SetParameters(NIterations, NBinsX, NBinsY, HTNBinsX, HTNBinsY, HTCut, HitCut, MaxDistance, MinDistance, MaxRadius, MinRadius, NParts, SqLen); CbmRichFindRings* richFindRings = new CbmRichFindRings(); richFindRings->UseFinder(richFinder); fRun->AddTask(richFindRings); //-------------------------------------------------------------------------- //-------------------- RICH Ring Fitting ----------------------------------- Double_t iRingCorr = 1.; // correction done (default), choose 0 if not CbmRichRingFitter* richFitter = new CbmRichRingFitterTAU(iVerbose, iRingCorr, 2); CbmRichFitRings* fitRings = new CbmRichFitRings("","",richFitter); fRun->AddTask(fitRings); //-------------------------------------------------------------------------- // ------------------- RICH Ring matching --------------------------------- CbmRichMatchRings* matchRings = new CbmRichMatchRings(iVerbose); fRun->AddTask(matchRings); // ------------------------------------------------------------------------- //--------------------- RICH ring-track assignment ------------------------ Double_t richDistance = 10.; // Max. dist. ring centre to track [cm] Int_t richNPoints = 5; // Minmum number of hits on ring CbmRichRingTrackAssign* richAssign = new CbmRichRingTrackAssignClosestD(richDistance, richNPoints, iVerbose); CbmRichAssignTrack* assignTrack = new CbmRichAssignTrack(); assignTrack->UseAssign(richAssign); fRun->AddTask(assignTrack); // ------------------------------------------------------------------------ //--------------------- RICH ring selection ------------------------------- CbmRichRingSelect2DCuts *ringSelect2D = new CbmRichRingSelect2DCuts(iVerbose); CbmRichSelectRings* richSelectRings2D = new CbmRichSelectRings(); richSelectRings2D->UseSelect(ringSelect2D); fRun->AddTask(richSelectRings2D); CbmRichRingSelectNeuralNet *ringSelectNN = new CbmRichRingSelectNeuralNet(iVerbose, "NeuralNet_RingSelection_Weights.txt"); CbmRichSelectRings* richSelectRingsNN = new CbmRichSelectRings(); richSelectRingsNN->UseSelect(ringSelectNN); fRun->AddTask(richSelectRingsNN); // ------------------------------------------------------------------------ // === End of RICH local reconstruction === // ========================================================================= // ========================================================================= // === ECAL reconstruction === // ========================================================================= // ----- ECAL hit producer ---------------------------------------------- CbmEcalHitProducerFastMC* ecalHitProd = new CbmEcalHitProducerFastMC("ECAL Hitproducer"); fRun->AddTask(ecalHitProd); // ------------------------------------------------------------------------- // === End of ECAL reconstruction === // ========================================================================= // ----- Parameter database -------------------------------------------- TString stsDigi = gSystem->Getenv("VMCWORKDIR"); stsDigi += "/parameters/sts/"; stsDigi += stsDigiFile; CbmRuntimeDb* rtdb = fRun->GetRuntimeDb(); CbmParRootFileIo* parIo1 = new CbmParRootFileIo(); CbmParAsciiFileIo* parIo2 = new CbmParAsciiFileIo(); parIo1->open(parFile.Data()); parIo2->open(stsDigi.Data(),"in"); rtdb->setFirstInput(parIo1); rtdb->setSecondInput(parIo2); rtdb->setOutput(parIo1); rtdb->saveOutput(); // ------------------------------------------------------------------------ // ----- Load transport geometry -------------------------------------- fRun->LoadGeometry(); // ------------------------------------------------------------------------ // ----- Intialise and run -------------------------------------------- fRun->Init(); fRun->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; // ------------------------------------------------------------------------ }