// -------------------------------------------------------------------------- // // Macro for reconstruction of simulated events // // All hitproducers (STS, RICH, TRD, TOF, ECAL), Fast MC for ECAL // STS track finding and fitting // TRD track finding and fitting // 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 // // -------------------------------------------------------------------------- { // ======================================================================== // Adjust this part according to your requirements // Verbosity level (0=quiet, 1=event level, 2=track level, 3=debug) Int_t iVerbose = 1; // Input file (MC events) TString inFile = "test.mc.root"; // Number of events to process Int_t nEvents = 2; // Background file (for MAPS pileup) TString bgFile = "auau.25gev.mbias.mc.root"; // Parameter file TString parFile = "params.root"; // Output file TString outFile = "test.reco.root"; // Version of STS digitisation scheme TString stsDigiVersion = "v05a"; // Specify MAPS digitisation parameters Int_t iMapsF = 1; // First MAPS station Int_t iMapsL = 2; // Last MAPS station Double_t mapsSigmaX = 0.0005; // resolution in x [cm] Double_t mapsSigmaY = 0.0005; // resolution in y [cm] Double_t mapsEff = 0.99; // single hit efficiency Double_t mapsFakeRate = 0.03; // fake hit rate Double_t mapsMergeDist = 0.00; // distance for hit merging [cm] Int_t mapsPileUp = 1; // number of pileup events // Specify Hybrid digitisation parameters Int_t iHybF = 3; // First hybrid station Int_t iHybL = 4; // Last hybrid station Double_t hybLx = 0.0050; // Pixel x size [cm] Double_t hybLy = 0.0050; // Pixel y size [cm] // Specify Strip digitisation parameters Int_t iStripF = 5; // First strip station Int_t iStripL = 8; // Last strip station // Specify parameters for RICH track extrapolation Int_t richNSts = 6; // minimum number of STS hits for extrapolation Double_t richZPos = 300.; // z position for extrapolation [cm] // Specify parameters for RICH projektion to PM plane Int_t richZFlag = 1; // Projetion from IM plane (default) // Specify parameters for RICH hit producer Double_t richPmtRad = 0.4; // PMT radius [cm] Double_t richPmtDist = 0.; // Distance between PMTs [cm] Int_t richDetType = 4; // Detector type Hammamatsu Int_t richNoise = 220; // Number of noise points per event // Specify whether ring radius correction to be done after ring fitting Double_t iRingCorr = 1.; // correction done (default), choose 0 if not // Specify TRD digitisation parameters 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, default = 100 Float_t trdDFoils = 0.001; // thickness of 1 foil [cm], default = 0.001 Float_t trdDGap = 0.008; // thickness of gap between foils [cm], default = 0.008 // Specify parameters for TOF digitisation Double_t tofSigmaT = 80.; Double_t tofSigmaXY = 0.5; // Specify parameters for RICH ring-track assignment Double_t richDistance = 1.; // Maximum distance between ring centre and track [cm] Int_t richNPoints = 5; // Minmum number of hits on ring // ---- Load libraries ------------------------------------------------- gROOT->LoadMacro("$VMCWORKDIR/gconfig/basiclibs.C"); basiclibs(); gSystem->Load("libGeoBase"); gSystem->Load("libParBase"); gSystem->Load("libBase"); gSystem->Load("libField"); gSystem->Load("libGen"); gSystem->Load("libPassive"); gSystem->Load("libSts"); gSystem->Load("libRich"); gSystem->Load("libTrd"); gSystem->Load("libTof"); gSystem->Load("libEcal"); gSystem->Load("libGlobal"); gSystem->Load("libKF"); gSystem->Load("libL1"); // ------------------------------------------------------------------------ // --- Now choose concrete engines for the different tasks ------------- CbmKF* kalman= new CbmKF(); CbmL1 *l1 = new CbmL1(); CbmStsTrackFinder* stsTrackFinder = new CbmL1StsTrackFinder(); CbmStsTrackFitter* stsTrackFitter = new CbmStsKFTrackFitter(); CbmTrdTrackFinder* trdTrackFinder = new CbmL1TrdTrackFinderSts(); ((CbmL1TrdTrackFinderSts*)trdTrackFinder)->SetVerbose(iVerbose); // ((CbmL1TrdTrackFinderSts*)trdTrackFinder)->SetPid(211); // ((CbmL1TrdTrackFinderSts*)trdTrackFinder)->SetChi2Cut(35); CbmTrdTrackFitter* trdTrackFitter = new CbmTrdTrackFitterKF(); ((CbmTrdTrackFitterKF*)trdTrackFitter)->SetVerbose(iVerbose); ((CbmTrdTrackFitterKF*)trdTrackFitter)->SetPid(211); CbmTrackMerger* trackMerger = new CbmTrackMergerIdeal(); CbmPrimaryVertexFinder* pvFinder = new CbmPVFinderIdeal(); CbmRichTrackExtrapolation* richExtra = new CbmRichTrackExtrapolationIdeal(richNSts, iVerbose); CbmRichRingFinder* richFinder = new CbmRichRingFinderIdeal(iVerbose); CbmRichRingFitter* richFitter = new CbmRichRingFitterCircle(iVerbose,iRingCorr); CbmRichRingTrackAssign* richAssign = new CbmRichRingTrackAssignClosestD(richDistance, richNPoints, iVerbose); // ------------------------------------------------------------------------ // In general, the following parts need not be touched // ======================================================================== // ----- Timer -------------------------------------------------------- TStopwatch timer; timer.Start(); // ------------------------------------------------------------------------ // ----- Reconstruction run ------------------------------------------- CbmRunAna *fRun= new CbmRunAna(); fRun->SetInputFile(inFile); fRun->SetOutputFile(outFile); // ------------------------------------------------------------------------ // ----- Parameter database -------------------------------------------- CbmRuntimeDb* rtdb = fRun->GetRuntimeDb(); CbmParRootFileIo* parInput1 = new CbmParRootFileIo(); parInput1->open(parFile.Data()); CbmParAsciiFileIo* parInput2 = new CbmParAsciiFileIo(); TString stsDigiFile = gSystem->Getenv("VMCWORKDIR"); stsDigiFile += "/parameters/sts/sts_digi_new_standard.par"; parInput2->open(stsDigiFile.Data(),"in"); rtdb->setFirstInput(parInput1); rtdb->setSecondInput(parInput2); fRun->LoadGeometry(); // ------------------------------------------------------------------------ // ========================================================================= // ====== Hit Producers ====== // ========================================================================= // ----- STS hit producers -------------------------------------------- CbmStsMapsHitProducer* mapsHitProd = new CbmStsMapsHitProducer(iMapsF, iMapsL, mapsSigmaX, mapsSigmaY, mapsEff, mapsFakeRate, mapsMergeDist, mapsPileUp, stsDigiVersion); mapsHitProd->SetPileupFileName(bgFile); CbmStsHybridHitProducer* hybHitprod = new CbmStsHybridHitProducer(iHybF, iHybL, hybLx, hybLy, stsDigiVersion); CbmStsStripHitProducer* stripHitProd = new CbmStsStripHitProducer(iStripF, iStripL, stsDigiVersion); fRun->AddTask(mapsHitProd); fRun->AddTask(hybHitprod); fRun->AddTask(stripHitProd); // ------------------------------------------------------------------------ // ---------------------RICH Hit Producer --------------------------------- CbmRichHitProducer* richHitProd = new CbmRichHitProducer(richPmtRad, richPmtDist, richDetType, richNoise, iVerbose); fRun->AddTask(richHitProd); //------------------------------------------------------------------------- // ----- TRD hit producer --------------------------------------------- CbmTrdHitProducer* trdHitProd = new CbmTrdHitProducer("TRD Hitproducer", "TRD task"); trdHitProd->SetPar(trdNFoils, trdDFoils, trdDGap); trdHitProd->SetSigmaX(trdSigmaX); trdHitProd->SetSigmaY(trdSigmaY1, trdSigmaY2, trdSigmaY3); fRun->AddTask(trdHitProd); // ------------------------------------------------------------------------ /* // ----- TOF hit producer --------------------------------------------- CbmTofHitProducer* tofHitProd = new CbmTofHitProducer("TOF HitProducer", iVerbose); tofHitProd->SetSigmaT(tofSigmaT); tofHitProd->SetSigmaXY(tofSigmaXY); fRun->AddTask(tofHitProd); // ------------------------------------------------------------------------ */ // ----- ECAL hit producer --------------------------------------------- CbmEcalHitProducerFastMC* ecalHitProd = new CbmEcalHitProducerFastMC("ECAL Hitproducer"); fRun->AddTask(ecalHitProd); // ------------------------------------------------------------------------ // ===== End of HitProducers ===== // ========================================================================= // ========================================================================= // ===== STS local tracking ====== // ========================================================================= // ----- STS track finding -------------------------------------------- CbmStsFindTracks* stsFindTracks = new CbmStsFindTracks("Track Finder", "CbmTask", stsTrackFinder); fRun->AddTask(kalman); fRun->AddTask(l1); fRun->AddTask(stsFindTracks); // ------------------------------------------------------------------------ // ----- STS track fitting -------------------------------------------- CbmStsFitTracks* stsFitTracks = new CbmStsFitTracks("STS Track Fitter", "CbmTask", stsTrackFitter); fRun->AddTask(stsFitTracks); // ------------------------------------------------------------------------ // ----- STS track matching ------------------------------------------- CbmStsMatchTracks* stsMatchTracks = new CbmStsMatchTracks("Match tracks", "STS", iVerbose); fRun->AddTask(stsMatchTracks); // ------------------------------------------------------------------------ // ===== End of STS local tracking ===== // ========================================================================= // ========================================================================= // ===== TRD local tracking ====== // ========================================================================= // ----- TRD track finding -------------------------------------------- CbmTrdFindTracks* trdFindTracks = new CbmTrdFindTracks("TRD Track Finder"); trdFindTracks->UseFinder(trdTrackFinder); fRun->AddTask(trdFindTracks); // ------------------------------------------------------------------------ // ----- TRD track fitting -------------------------------------------- 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 tracking ===== // ========================================================================= // ========================================================================= // ===== Global tracking ====== // ========================================================================= // ----- Global tracking ---------------------------------------------- CbmFindGlobalTracks* globalFinder = new CbmFindGlobalTracks(trackMerger, NULL, NULL, iVerbose); fRun->AddTask(globalFinder); // ------------------------------------------------------------------------ // ----- Primary vertex finding --------------------------------------- CbmFindPrimaryVertex* pvFindTask = new CbmFindPrimaryVertex(pvFinder); fRun->AddTask(pvFindTask); // ------------------------------------------------------------------------ // ===== End of global tracking ===== // ========================================================================= // ========================================================================= // ===== RICH pattern recognition ====== // ========================================================================= //----------------------RICH Track Extrapolation -------------------------- CbmRichExtrapolateTracks* richExtrapolate = new CbmRichExtrapolateTracks(); richExtrapolate->UseExtrapolation(richExtra,richZPos); fRun->AddTask(richExtrapolate); //-------------------------------------------------------------------------- //--------------------- Rich Track Projection to photodetector ------------- CbmRichProjectionProducer* richProj = new CbmRichProjectionProducer(iVerbose, richZFlag); fRun->AddTask(richProj); //-------------------------------------------------------------------------- //----------------------RICH Ring Finding---------------------------------- CbmRichFindRings* richFindRings = new CbmRichFindRings(); richFindRings->UseFinder(richFinder); fRun->AddTask(richFindRings); //------------------------------------------------------------------------- //-------------------- RICH Ring Fitting ---------------------------------- CbmRichFitRings* fitRings = new CbmRichFitRings("","",richFitter); fRun->AddTask(fitRings); //------------------------------------------------------------------------- // ------------------- RICH Ring matching -------------------------------- CbmRichMatchRings* matchRings = new CbmRichMatchRings(iVerbose); fRun->AddTask(matchRings); // ------------------------------------------------------------------------ //--------------------- RICH ring-track assignment ------------------------ CbmRichAssignTrack* assignTrack = new CbmRichAssignTrack(); assignTrack->UseAssign(richAssign); fRun->AddTask(assignTrack); // ------------------------------------------------------------------------ // ===== End of RICH pattern recognition ===== // ========================================================================= // ----- 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; // ------------------------------------------------------------------------ }