outdated/dc1/stt/runreco.C

Go to the documentation of this file.

{
// ========================================================================
// Verbosity level (0=quiet, 1=event level, 2=track level, 3=debug)
Int_t iVerbose = 0;

// Input file (MC events)
TString inFile = "points_sttmvdemc.root";

// Parameter file
TString parFile = "testparams.root";

// Output file
TString outFile = "tracks_sttmvdemc.root";

// Number of events to process
Int_t nEvents = 0;
 
// ----  Load libraries   -------------------------------------------------
gROOT->LoadMacro("$VMCWORKDIR/gconfig/basiclibs.C");
basiclibs();
gSystem->Load("libGeoBase");
gSystem->Load("libParBase");
gSystem->Load("libBase");
gSystem->Load("libPndData");
gSystem->Load("libField");
gSystem->Load("libPassive");
gSystem->Load("libGen");  
gSystem->Load("libStt"); 
gSystem->Load("libEmc"); 
gSystem->Load("libgenfit");
gSystem->Load("libtpc"); 
gSystem->Load("libtpcreco");
gSystem->Load("libtrackrep");
gSystem->Load("librecotasks");
gSystem->Load("libMvd");
gSystem->Load("libMvdReco");
gSystem->Load("libLHETrack");
// ------------------------------------------------------------------------

// ---  Now choose concrete engines for the different tasks   -------------
// ------------------------------------------------------------------------

// In general, the following parts need not be touched
// ========================================================================

// -----   Timer   --------------------------------------------------------
TStopwatch timer;
timer.Start();
// ------------------------------------------------------------------------

// -----   Digitization run   -------------------------------------------
FairRunAna *fRun= new FairRunAna();
fRun->SetInputFile(inFile);
fRun->SetOutputFile(outFile);
// ------------------------------------------------------------------------

// -----  Parameter database   --------------------------------------------
TString allDigiFile = "./all.par";

FairRuntimeDb* rtdb = fRun->GetRuntimeDb();
FairParRootFileIo* parInput1 = new FairParRootFileIo();
parInput1->open(parFile.Data());
        
FairParAsciiFileIo* parIo1 = new FairParAsciiFileIo();
parIo1->open(allDigiFile.Data(),"in");
        
rtdb->setFirstInput(parInput1);
rtdb->setSecondInput(parIo1);
fRun->LoadGeometry();
// ------------------------------------------------------------------------

// -----   STT analysis tasks   --------------------------------------------
//PndSttHitProducerIdeal* sttHitProducer = new PndSttHitProducerIdeal();
PndSttHitProducerRealFast* sttHitProducer = new PndSttHitProducerRealFast();
fRun->AddTask(sttHitProducer);

// trackfinding ....
PndSttTrackFinderIdeal* sttTrackFinder = new PndSttTrackFinderIdeal(iVerbose);
PndSttFindTracks* sttFindTracks = new PndSttFindTracks("Track Finder", "FairTask", sttTrackFinder, iVerbose);
sttFindTracks->AddHitCollectionName("STTHit", "STTPoint");
fRun->AddTask(sttFindTracks);

// trackmatching ....
PndSttMatchTracks* sttTrackMatcher = new PndSttMatchTracks("Match tracks", "STT", iVerbose);
sttTrackMatcher->AddHitCollectionName("STTHit", "STTPoint");
fRun->AddTask(sttTrackMatcher);  

// trackfitting ....
PndSttTrackFitter* sttTrackFitter = new PndSttHelixTrackFitter(0);
PndSttFitTracks* sttFitTracks = new PndSttFitTracks("STT Track Fitter", "FairTask", sttTrackFitter); 
sttFitTracks->AddHitCollectionName("STTHit");
fRun->AddTask(sttFitTracks);

// -----   MDV digi producers   --------------------------------- 
// DIGI
double   topPitch=0.015,//cm
  botPitch=0.015,//cm
  orient=TMath::Pi()*(0.5),
  skew=TMath::Pi()*(0.5);
TVector2 topAnchor(0.,0.);
TVector2 botAnchor(0.,0.);
int      topFE = 10,
  botFE = 4,
  nrFEChannels = 128;
double   threshold=3000., noise=1000.;

PndMvdStripHitProducer* mvdHitProd
= new PndMvdStripHitProducer(topPitch, botPitch,
                             orient, skew,
                             topAnchor,botAnchor,
                             topFE, botFE, nrFEChannels,
                             threshold, noise);
mvdHitProd->SetVerbose(iVerbose);
fRun->AddTask(mvdHitProd);

Double_t  lx=0.01, ly=0.01, threshold=600, noise=200;
PndMvdHybridHitProducer* mvdPixProd = new PndMvdHybridHitProducer(lx,ly,threshold,noise);
mvdPixProd->SetVerbose(iVerbose);
fRun->AddTask(mvdPixProd);
 
// CLUST
// Cluster finding for strip detectors
Double_t noise = 1000.; 
Double_t chargecut = 3. * noise;
PndMvdStripClusterTask* mvdmccls = new PndMvdStripClusterTask(chargecut,inFile);
mvdmccls->SetVerbose(iVerbose);
fRun->AddTask(mvdmccls);
  
// Cluster finder for pixel detectors
PndMvdPixelClusterTask* mvdClusterizer = new PndMvdPixelClusterTask(1.8,76,84, inFile);
mvdClusterizer->SetVerbose(iVerbose);
fRun->AddTask(mvdClusterizer);
 
// TRACKFINDER
//PndMvdIdealTrackingTask* mvdmctrk = new PndMvdIdealTrackingTask();
//mvdmctrk->SetVerbose(iVerbose);
//fRun->AddTask(mvdmctrk);

// -----   EMC hit producers   ---------------------------------
PndEmcHitProducer* emcHitProd = new PndEmcHitProducer();
fRun->AddTask(emcHitProd); // hit production 

PndEmcMakeDigi* emcMakeDigi=new PndEmcMakeDigi();
//fRun->AddTask(emcMakeDigi); // fast digitization

PndEmcHitsToWaveform* emcHitsToWaveform= new PndEmcHitsToWaveform(iVerbose);
PndEmcWaveformToDigi* emcWaveformToDigi=new PndEmcWaveformToDigi(iVerbose);
fRun->AddTask(emcHitsToWaveform);  // full digitization
fRun->AddTask(emcWaveformToDigi);  // full digitization

PndEmcMakeCluster* emcMakeCluster= new PndEmcMakeCluster(iVerbose);
fRun->AddTask(emcMakeCluster);

PndEmcHdrFiller* emcHdrFiller = new PndEmcHdrFiller();
fRun->AddTask(emcHdrFiller); // ECM header

PndEmcMakeBump* emcMakeBump= new PndEmcMakeBump();
fRun->AddTask(emcMakeBump);

PndEmcMakeRecoHit* emcMakeRecoHit= new PndEmcMakeRecoHit();
fRun->AddTask(emcMakeRecoHit);

// -----   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;
// ------------------------------------------------------------------------


}