/* $Id: CbmRichOnlyReconstruction.C,v 1.5 2006/09/15 12:50:52 turany Exp $ */ /* History of CVS commits: * * $Log: CbmRichOnlyReconstruction.C,v $ * Revision 1.5 2006/09/15 12:50:52 turany * load the base before the field * * Revision 1.4 2006/06/22 07:18:07 hoehne * libraries adopted to recent changes * * Revision 1.3 2006/04/10 15:21:02 hoehne * libraries had to be reorderd after recent changes in CVS * * Revision 1.2 2006/04/10 09:52:32 hoehne * more ring finders added with example how to run * * Revision 1.1 2006/04/07 16:09:09 hoehne * initial version: RICH only reconstruction, tracking to be required from a separate root file * * * */ // -------------------------------------------------------------------------- // // Macro for RICH reconstruction // to be run after STS reconstruction // input: file with MC detector simulation (+parfile) // file with STS track reconstruction // output: file with RICH reconstruction // // note: global track cannot be filled with RICH information this way! // // 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 // // // C. Hoehne // // -------------------------------------------------------------------------- { // ======================================================================== // 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/auau.25gev.centr.mc_notrd_test2.root"; //TString inFile ="/d/cbm02/hoehne/data/urqmd/auau.25gev.centr.mc_notrd.root"; //TString inFile ="/d/cbm02/slebedev/rich/auau.25gev.centr.mc.0001.root"; //TString inFile2 ="/d/cbm02/slebedev/rich/auau.25gev.centr.reco.0001.root"; TString inFile ="../run/test.mc.root"; TString inFile2 ="../rich/test.reco.rich.root"; // Number of events to process Int_t nEvents = 10; // Parameter file //TString parFile = "../run/parfiles/params_N2_activeMagnet.root"; //TString parFile = "../run/parfiles/params_N2_activeMagnet_noPMTglass_test2.root"; //TString parFile = "/d/cbm02/slebedev/rich/parfiles/params.0001.root"; TString parFile = "../run/params.root"; // Output file TString outFile = "test.reco.rich2.root"; //TString outFile = "test.root"; //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=CsI, 3=Hamamatsu) 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->AddFriend(inFile2); 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(); // ------------------------------------------------------------------------ // Kalman Filter Instance (needed for running the ExtrapolationKF) CbmKF *KF = new CbmKF(); fRun->AddTask(KF); //----------------------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 CbmRichTrackExtrapolationKF* extrapolation = new CbmRichTrackExtrapolationKF(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 // Ideal Ring Finder // Hough Transform // Elastic Net // track based finder CbmRichRingFinderIdeal* ringFinder = new CbmRichRingFinderIdeal(iVerbose); /* CbmRichRingFinderHough* ringFinder = new CbmRichRingFinderHough(iVerbose); //describtion of these parameters see in CbmRichRingFinderHough.cxx file int NIterations = 2; int NBinsX[] = {250,250,150};//550,550,350 int NBinsY[] = {100,100,60};//200,200,450 int HTNBinsX[] = {150,150,75};//150,150,100 int HTNBinsY[] = {60,60,30};//200,200,120 int HTCut[] = {200,100,50};//500,60,5 int HitCut[] = {20,15,8};//50,15,3 double MaxDistance[] = {12.0,12.0,12.0}; double MinDistance[] = {1.0,1.0,1.0}; double MinRadius[] = {2.0,2.0,2.0}; double MaxRadius[] = {6.0,6.0,6.0}; int NParts[] = {2,1,1}; int SqLen = 2; ringFinder.SetParameters(NIterations,NBinsX,NBinsY,HTNBinsX,HTNBinsY,HTCut,HitCut, MaxDistance,MinDistance,MaxRadius,MinRadius,NParts,SqLen); */ // CbmL1RichENNRingFinder * ringFinder = new CbmL1RichENNRingFinder(); /* CbmRichRingFinderTrack * ringFinder = new CbmRichRingFinderTrack ( 1 ); ringFinder->SetParameters(20,3.0,6.5,12); */ 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,2); CbmRichRingTrackAssignIdeal *Tassign = new CbmRichRingTrackAssignIdeal(distance,npoints,2); CbmRichAssignTrack* AssignTrack = new CbmRichAssignTrack(); AssignTrack->UseAssign(Tassign); fRun->AddTask(AssignTrack); //---------------------RICH quality checking of ring finder ------------------------------ /* CbmRichRingQa *RingQa = new CbmRichRingQa("CbmRichRingQa","CbmRichRingQa T",iVerbose); fRun->AddTask(RingQa); */ // ----- Intialise and run -------------------------------------------- fRun->Init(); fRun->Run(0,nEvents); // ------------------------------------------------------------------------ // RingQa->Finish(); // ----- 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; // ------------------------------------------------------------------------ }