// -------------------------------------------------------------------------- // // Macro for STS + TRD reconstruction of simulated events // full RichReconstruction added // // STS and TRD Hitproducers // L1 STS tracking routines // Ideal track finding in TRD // Matching of reconstructed and MC tracks in STS and TRD // Global tracking (STS+TRD), ideal // Ideal (dummy) primary vertex finder // // RICH: // track extrapolation to z-plane in Rich detector // track projection to photodetector plane // Rich Hit Producer // Rich Ring Finder // Rich Ring Fitter // Matching of reconstructed and MC rings // ring-track assignement // // // V. Friese 28/11/2005 // C. Hoehne 26/01/2006 RICH reconstruction added // // -------------------------------------------------------------------------- { // ======================================================================== // 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 = "../run/auau.25gev.centr.mc_100notrd.root"; TString inFile = "../run/test.mc.root"; //TString inFile ="/d/cbm02/hoehne/data/urqmd/auau.25gev.centr.mc_notrd.root"; // Number of events to process Int_t nEvents = 1; // Background file (for MAPS pileup) TString bgFile = "../run/test.mc.root"; // Parameter file //TString parFile = "../run/parfiles/params_N2_activeMagnet.root"; TString parFile = "../run/params.root"; // Output file TString outFile = "test.reco.rich.root"; //TString outFile = "test.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 = 0; // number of pileup events // Specify Hybrid digitisation parameters Int_t iHybF = 3; // First hybrid station Int_t iHybL = 5; // 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 = 6; // First strip station Int_t iStripL = 7; // Last strip station // 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 //RICH //------------------------------------------------------------------------- // ----------------- RICH Track Extrapolation ----------------------------- // ----- parameters for track selection Int_t Nsts = 6; // minimum number of STS hits // ----- parameters for track extrapolation Double_t zPos = 300.; // z-Position for extrapolation [cm] // Specify parameter for RICH projection to photodetector plane Int_t zflag = 1; // zflag = 1 projection from imaginary z-plane (default) // ATTENTION: never use CbmRichTrackExtrapolationMirrorIdeal in this case // zflag = 2 projection from point array in mirror // ATTENTION: in this case CbmRichTrackExtrapolationMirrorIdeal has to be run!! // ---------------------RICH Hit Producer --------------------------------- Double_t pmt_rad = 0.3; // PMT radius Double_t pmt_dist = 0.; // distance between PMTs Int_t det_type = 1; // detector type (choose: 1=Protvino, 2=Hamamatsu, 3=CsI) Int_t noise = 220; // number of noise points to be added // (note: excluding geom. eff. (~0.9 for r=0.3 and d=0)) //---------------------RICH ring-track assignement --------------------------------- // -----parameters ------------------------------------------------------------- Double_t distance = 1.; // maximum distnace between ring center and track [cm] Int_t npoints = 5; // minimum number of hits/ring // In general, the following parts need not be touched // ======================================================================== // ----- 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("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"); // ------------------------------------------------------------------------ // ----- 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.par"; parInput2->open(stsDigiFile.Data(),"in"); rtdb->setFirstInput(parInput1); rtdb->setSecondInput(parInput2); fRun->LoadGeometry(); // ------------------------------------------------------------------------ // ----- 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); // ------------------------------------------------------------------------ // ----- 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); // ------------------------------------------------------------------------ // ----- STS track finding -------------------------------------------- CbmKF* kalman= new CbmKF(); fRun->AddTask(kalman); CbmL1 *l1 = new CbmL1(); fRun->AddTask(l1); CbmStsTrackFinder* stsTrackFinder = new CbmL1StsTrackFinder(); CbmStsFindTracks* stsFindTracks = new CbmStsFindTracks("Track Finder", "CbmTask", stsTrackFinder); fRun->AddTask(stsFindTracks); // ------------------------------------------------------------------------ // ----- STS track matching ------------------------------------------- CbmStsMatchTracks* stsMatchTracks = new CbmStsMatchTracks("Match tracks", "STS", iVerbose); fRun->AddTask(stsMatchTracks); // ------------------------------------------------------------------------ // ----- STS track fitting -------------------------------------------- CbmStsTrackFitter* stsTrackFitter = new CbmStsKFTrackFitter(); CbmStsFitTracks* stsFitTracks = new CbmStsFitTracks("STS Track Fitter", "CbmTask", stsTrackFitter); fRun->AddTask(stsFitTracks); // ------------------------------------------------------------------------ /* // ----- TRD track finding -------------------------------------------- CbmL1TrdTrackFinderSts* trdTrackFinder = new CbmL1TrdTrackFinderSts(); trdTrackFinder->SetVerbose(iVerbose); trdTrackFinder->SetPid(211); trdTrackFinder->SetChi2Cut(35); CbmTrdFindTracks* trdFindTracks = new CbmTrdFindTracks("TRD Track Finder"); trdFindTracks->UseFinder(trdTrackFinder); fRun->AddTask(trdFindTracks); // ------------------------------------------------------------------------ // ----- TRD track matching ------------------------------------------- CbmTrdMatchTracks* trdMatchTracks = new CbmTrdMatchTracks(iVerbose); fRun->AddTask(trdMatchTracks); // ------------------------------------------------------------------------ // ----- TRD track fitting -------------------------------------------- CbmL1TrdTrackFitter* trdTrackFitter = new CbmL1TrdTrackFitter(); trdTrackFitter->SetVerbose(iVerbose); trdTrackFitter->SetPid(211); CbmTrdFitTracks* trdFitTracks = new CbmTrdFitTracks("TRD track fitter", "TRD", trdTrackFitter); // ------------------------------------------------------------------------ */ // ----- Global tracking ---------------------------------------------- CbmTrackMerger* trackMerger = new CbmTrackMergerIdeal(); CbmFindGlobalTracks* globalFinder = new CbmFindGlobalTracks(trackMerger, NULL, NULL, iVerbose); fRun->AddTask(globalFinder); // ------------------------------------------------------------------------ // ----- Primary vertex finding --------------------------------------- CbmPrimaryVertexFinder* pvFinder = new CbmPVFinderIdeal(); CbmFindPrimaryVertex* pvFindTask = new CbmFindPrimaryVertex(pvFinder); fRun->AddTask(pvFindTask); // ------------------------------------------------------------------------ //----------------------RICH Track Extrapolation --------------------------------- // This is the place to choose a particular track extrapolation - available are: // CbmRichTrackExtrapolationIdeal: select MCPoints for extrapolation (RichImPlanePoint) // CbmRichTrackExtrapolationKF : extrapolation of STS tracks with Kalman Fitter CbmRichTrackExtrapolationIdeal* extrapolation = new CbmRichTrackExtrapolationIdeal (Nsts,iVerbose); CbmRichExtrapolateTracks* richExtrapolate = new CbmRichExtrapolateTracks(); richExtrapolate->UseExtrapolation(extrapolation,zPos); fRun->AddTask(richExtrapolate); //-------------------------------------------------------------------------- //--------------------- Rich Track Projection to photodetector ------------- CbmRichProjectionProducer *proj= new CbmRichProjectionProducer(iVerbose,zflag); fRun->AddTask(proj); // ---------------------RICH Hit Producer --------------------------------- CbmRichHitProducer *hp= new CbmRichHitProducer(pmt_rad,pmt_dist,det_type,noise,iVerbose); fRun->AddTask(hp); //------------------------------------------------------------------------- //----------------------RICH Ring Finding---------------------------------- // This is the place to choose a particular finder CbmRichRingFinderIdeal* ringFinder = new CbmRichRingFinderIdeal(iVerbose); CbmRichFindRings* richFindRings = new CbmRichFindRings(); richFindRings->UseFinder(ringFinder); fRun->AddTask(richFindRings); //-------------------------------------------------------------------------- //--------------------RICH Ring Fitting------------------------------------- // This is the place to choose a particular fitting method - available are: // CbmRichRingFitterCircle // CbmRichRingFitterCOP // CbmRichRingFitterRobustCOP // CbmRichRingFitterTAU CbmRichRingFitterCircle* ringFitter = new CbmRichRingFitterCircle(); CbmRichFitRings* fitRings = new CbmRichFitRings("","",ringFitter); fRun->AddTask(fitRings); //-------------------------------------------------------------------------- // ------------------- RICH Ring matching ---------------------------------- CbmRichMatchRings* matchRings = new CbmRichMatchRings(iVerbose); fRun->AddTask(matchRings); // ------------------------------------------------------------------------ //---------------------RICH ring-track assignement --------------------------------- CbmRichRingTrackAssignClosestD *Tassign = new CbmRichRingTrackAssignClosestD(distance,npoints,iVerbose); CbmRichAssignTrack* AssignTrack = new CbmRichAssignTrack(); AssignTrack->UseAssign(Tassign); fRun->AddTask(AssignTrack); // ----- 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; // ------------------------------------------------------------------------ }