void runReco_vw(unsigned int run,
		TString outpath,
		TString outFileName = "NULL",
		bool clcorr=false,
		bool reclust=false,
		unsigned int nEvents = 0,
		unsigned int nEvStart = 0,
		unsigned int bfield=1,
		unsigned int smoothing = 1,
	        int outnum=-1, 
		unsigned int gaincorr=1,
		double fac = 4.0,
		double stepoff = 1.0,
		double fC=14,
		double fG=-1)
{
// ========================================================================
// Verbosity level (0=quiet, 1=event level, 2=track level, 3=debug)
  Int_t iVerbose = 0;

  TStopwatch timer;
  timer.Start();
  
  if (outFileName != "NULL")
    {
      TString outfilename = outFileName;
      outfilename += ".root";
      
      TString filename="CosmicData/";
      filename+="runC_";
      filename+=run;
      filename+=".lmd_decoded.root";
    }

  else
    {
      
      TString filename="decoded/";
      filename+="runC_";
      filename+=run;
      filename+=".lmd_decoded.root";
    }
  TString basedir = gSystem->Getenv("VMCWORKDIR");
  FairRunAna* fRun = new FairRunAna();

  TString jobdir = outpath; 
  std::string jobname(filename.Data());
  int last = jobname.rfind("/");
  if(last>0)
    jobname = jobname.substr(last+1,jobname.size()+1);  

  //create output file -----------------------------------------------  
  TString outName(jobname);
  if(outName.Contains("repaired."))
    outName.ReplaceAll(".lmd_decoded_repaired","");
  if(outName.Contains("decoded."))
    outName.ReplaceAll(".lmd_decoded","");
  if(smoothing)
    outName.ReplaceAll(".root", "_smoothed.root");
  if(clcorr)
    outName.ReplaceAll(".root","_clcorr.root");
  if (reclust)
    outName.ReplaceAll(".root","_recl.root");
  if(outnum>-1)
    outName.ReplaceAll(".root",Form("_%03i.root",outnum));


  TString outFile = outpath+"/";
  if (outFileName != "NULL")
    {
      cout<<outFile<<outfilename<<"\n";
      outFile += outfilename;
    }
  else
    {
      outFile += outName;
    }

  outName.ReplaceAll(".root",".reco.root");
  cout<<"OUTFILE:"<<outFile<<"\n";
  fRun->SetOutputFile(outFile);
  //DONE create output file -----------------------------------------------  
 
  FairRuntimeDb* rtdb = fRun->GetRuntimeDb();
  FairParAsciiFileIo* parInput1 = new FairParAsciiFileIo();
  TString parFile = basedir;
  parFile+="/tpc/FOPI/par/tpc.run";
  parFile+=run;
  parFile +=".par";
  cout<<parFile<<endl;
  fstream fin;
  fin.open(parFile.Data(),ios::in);
  if( !fin.is_open() )
    {
      cout<<"parfile doesnt exist"<<endl;
      return;
    }
  fin.close();
  parInput1->open(parFile.Data(),"in"); 
  rtdb->setFirstInput(parInput1);

  TString parOutFile=outFile;
  parOutFile.ReplaceAll("reco.root","param.root");
  FairParRootFileIo* parOut=new FairParRootFileIo(kTRUE);
  parOut->open(parOutFile.Data());
  //rtdb->setOutput(parOut);

  TpcDigiPar* par = (TpcDigiPar*) rtdb->getContainer("TpcDigiPar");
  par->setInputVersion(fRun->GetRunId(),1);
  par->setChanged(kTRUE);
  rtdb->Print();  
  par->init();


  Double_t taroff = -65;
  Double_t targetpos = taroff+41.;
  // possible materials: copper (default), lead, carbon
  TString targetmat = TString("carbon");
  TString alignmentFileName = par->getAlignmentFile();
  TString geoFile;
  geoFile.Form("$FOPI2ROOT/fopigeometry/FopiGeom_s339%5.1f.root",targetpos);
  TFile* geotmp =new TFile(geoFile.Data(),"READ");
  if(geotmp->IsOpen())
    {
      std::cout<<"geofile exists already\n";
      geotmp->Close();
    }
  else
    {
      TString cmd;
      cmd.Form(".! root -b -q '$FOPI2ROOT/fopigeometry/FopiGeom.C(%f,\"%s\",\"%s\",\"%s\",kFALSE)'",targetpos,targetmat.Data(),alignmentFileName.Data(),geoFile.Data());
      gROOT->ProcessLine(cmd);
    }
  // add geometry file to fRun
  fRun->SetGeomFile(geoFile);
   

  FOPIField *fMagField = new FOPIField(0.616);
  fMagField->SetTargetOffset(taroff);
  fRun->SetField(fMagField);

  //initialize GF field
  GFFieldManager::getInstance()->init(new PndFieldAdaptor(fRun->GetField()));
  GFMaterialEffects::getInstance()->init(new GFTGeoMaterialInterface());

  /*
  //extract number of entries in external data tree
  cout<<"filename: "<<filename<<endl;
  TFile testFile(filename);
  TTree* testtree = testFile.Get("tpcEvent");
  unsigned int numEvents = testtree->GetEntries();
  std::cout<<"Found "<<numEvents<<" events in input data file"<<std::endl;
  if(nEvents==0 || nEvents>numEvents) nEvents=numEvents;
  testFile.Close();
  */
  //--------------------SET UP TASKS ------------------------------
  TpcEventCounter* evCount =  new TpcEventCounter();
  evCount->SetnEvts(nEvents-nEvStart);
  //evCount->SetStep(1);
  fRun->AddTask(evCount); 


  TpcDataReaderTask* read = new TpcDataReaderTask();
  read->SetPersistence();
  read->SetDatafile(filename);
  read->SetSampleBranchName("TpcSample");
  read->SetStartEvent(nEvStart);
  //read->SetCutSmallPad();
  //read->SetMinSamples(1000);
  fRun->AddTask(read);
  

  TpcPSATask* tpsa = new  TpcPSATask();
  tpsa->SetPersistence();
  //tpsa->SetSamplePersistence();
  tpsa->SetSampleBranchName("TpcSample"); // Input of PSA
  fRun->AddTask(tpsa);


  if (gaincorr==1)
    {
      TpcDigiAmpCorrectionTask* tpcCalib = new TpcDigiAmpCorrectionTask();
      fRun->AddTask(tpcCalib);
    }
  
  TpcClusterFinderTask* tpcCF = new TpcClusterFinderTask();
  if (clcorr)
    tpcCF->SetClusterBranchName("TpcCluster_raw");
  else if (reclust)
    tpcCF->SetClusterBranchName("TpcPreCluster");
  else
    tpcCF->SetClusterBranchName("TpcCluster");
  tpcCF->SetDigiPersistence(); // keep Digis copys in clusters
  tpcCF->SetPersistence(); // keep Clusters
  tpcCF->timeslice(6); //in samples
  tpcCF->SetErrorPars(-1,14);
  tpcCF->SetSingleDigiClusterAmpCut(0.);
  tpcCF->SetClusterAmpCut(0.);
  tpcCF->SetVerbose(iVerbose);
  tpcCF->SetSimpleClustering(); // use TpcClusterFinderSimple
  fRun->AddTask(tpcCF);

  if (clcorr)
    {
      cout<<"Setting up cluster correction task"<<endl;
      TpcClusterCorrectionTask* tpcCC = new TpcClusterCorrectionTask();
      tpcCC->SetClusterBranchName("TpcCluster_raw");
      if (reclust)
	tpcCC->SetClusterOutBranchName("TpcPreCluster");
      else
	tpcCC->SetClusterOutBranchName("TpcCluster");
      tpcCC->SetDevMap("datafiles/full_field_-148.93_-312.89_4s3_4.7__merged7.txt");
      fRun->AddTask(tpcCC);
    }
  
  cout<<"Setting up Rieman Task"<<endl;
  //find TpcRiemannTracks in the TPC alone
  TpcRiemannTrackingTask* tpcSPR1 = new TpcRiemannTrackingTask();
  if (reclust)
    tpcSPR1->SetClusterBranchName("TpcPreCluster");
  else
    tpcSPR1->SetClusterBranchName("TpcCluster");
  tpcSPR1->SetRiemannTrackBranchName("RiemannTrack");
  tpcSPR1->SetRiemannHitBranchName("RiemannHit");
  tpcSPR1->SetPersistence();
  tpcSPR1->SetSortingParameters(true,1,0);
  tpcSPR1->SetTrkFinderParameters(2.0,1.6,3,0.4);
  //tpcSPR1->SetTrkMergerParameters(16.,0.4,0.2,0.5)
  fRun->AddTask(tpcSPR1);

  if (reclust)
    {
      cout<<"Setting up TrackInitTask"<<endl;
      //build GFTracks from TpcRiemannTracks
      TpcTrackInitTask* trackInit1=new TpcTrackInitTask();
      trackInit1->SetPersistence();
      trackInit1->SetRiemannBranchName("RiemannTrack");
      trackInit1->SetTrackOutBranchName("TpcPreTrackPreFit");
      trackInit1->SetClusterBranchName("TpcPreCluster");
      trackInit1->SetRecoHitOutBranchName("TpcPreSPHit");
      trackInit1->SetCosmicSorting();
      //trackInit1->SetVerbose(iVerbose);
      trackInit1->SetPDG(13);
      //trackInit1->SetAlign(kFALSE);
      if (smoothing==1)
	trackInit1->SetSmoothing(true);
      fRun->AddTask(trackInit1);


      cout<<"Setting up KalmanTask 1"<<endl;
      KalmanTask* kalman1 =new KalmanTask();
      kalman1->SetTpcClusterBranchName("TpcPreSPHit");
      kalman1->SetTrackBranchName("TpcPreTrackPreFit");
      kalman1->SetOutBranchName("TrackPrePostFit");
      kalman1->SetPersistence();
      kalman1->SetNumIterations(1); // number of fitting iterations (back and forth)
      fRun->AddTask(kalman1);

      cout<<"Setting up Reclusterizer"<<endl;
      TpcReClusterizerTask* tpcReCl = new TpcReClusterizerTask();
      tpcReCl->SetTrackBranchName("TrackPrePostFit");
      tpcReCl->SetClusterBranchName("TpcPreCluster");
      tpcReCl->SetRecoHitBranchName("TpcPreSPHit");
      tpcReCl->SetReClusterBranchName("TpcCluster");
      //tpcReCl->SetAlign(kFALSE);
      //tpcReCl->SetStepSize(1);
      tpcReCl->SetSkipUnFitted();
      tpcReCl->SetMinClusterSize(3);
      tpcReCl->SetErrorPars(fG,fC);
      tpcReCl->SetVarStepSqrt();
      tpcReCl->SetStepOff_sqrt(stepoff);
      tpcReCl->SetFac(fac);
      tpcReCl->SetStepOff_0_sqrt(0);
      //tpcReCl->SetVerbose(iVerbose);
      //tpcReCl->SetTimeInfo();
      tpcReCl->SetPersistence();
      tpcReCl->SetAllRecoHitPersistence();
      tpcReCl->SetUseFirstDigiPos();
      //tpcReCl->SetUseChamberEdge();
      tpcReCl->SetUseCosmics();
      //tpcReCl->SetUseAllDigis();
      if(clcorr)
	tpcReCl->SetClusterCorr("datafiles/full_field_-148.93_-312.89_4s3_4.7__merged7.txt");
      tpcReCl->SetTrackInit();
      fRun->AddTask(tpcReCl);
    }
  else
    {
      //build GFTracks from TpcRiemannTracks
      TpcTrackInitTask* trackInit=new TpcTrackInitTask();
      trackInit->SetRiemannBranchName("RiemannTrack");
      trackInit->SetClusterBranchName("TpcCluster");
      trackInit->SetTrackOutBranchName("TpcTrackPreFit");
      trackInit->SetRecoHitOutBranchName("TpcSPHit");
      trackInit->SetPersistence();
      trackInit->SetCosmicSorting();
      trackInit->SetVerbose(iVerbose);
      //trackInit->SetMCPid(); // use ideal particle identification
      trackInit->SetPDG(13);
      //trackInit->useGeane(); // uses RKTrackrep and GeaneTrackrep
      if (smoothing==1)
	trackInit->SetSmoothing(true);
      fRun->AddTask(trackInit);
    }
  /*
  TpcAlignmentTask* align = new TpcAlignmentTask();
  align->SetRecoHitBranchName("TpcSPHit");
  fRun->AddTask(align);
  */

  KalmanTask* kalman =new KalmanTask();
  kalman->SetTpcClusterBranchName("TpcSPHit");
  kalman->SetTrackBranchName("TpcTrackPreFit");
  kalman->SetOutBranchName("TrackPostFit");
  kalman->SetPersistence();
  //kalman->SetClusterBranchName("TpcCluster_cut");
  kalman->SetNumIterations(3); // number of fitting iterations (back and forth)
  fRun->AddTask(kalman);

  cout<<"Setting up CosmicsTask"<<endl;
  TpcCosmicsTask* Cos = new TpcCosmicsTask();
  Cos->SetClusterBranchName("TpcCluster");
  Cos->SetPersistence();
  Cos->SetNumberOfTrackReps(1); // set to 2 if you use GeaneTrackrep (tpcSPR->useGeane();)
  Cos->SetDatafile(filename);
  if (smoothing==1)
    Cos->SetUnbiased();
  //Cos->SetVerbose();
  fRun->AddTask(Cos);

  cout<<"Setting up Track Cleaner"<<endl;
  TpcTrackCleaner* tpcTrCl = new TpcTrackCleaner();
  tpcTrCl->SetCosmics();
  tpcTrCl->SetThetaCut(40,140);
  tpcTrCl->SetClMeanCut(900);
  //tpcTrCl->SetPhiCut(70,110);
  //tpcTrCl->SetVerbose(iVerbose);
  tpcTrCl->SetKeepDelTracks();
  fRun->AddTask(tpcTrCl);  

  // -----   Intialise and run   --------------------------------------------
 
  
  cout<<"************start ini***************"<<endl;
  fRun->Init();
  //Initialize te Digimapper:
  TpcDigiPar* tpcpar = FairRun::Instance()->GetRuntimeDb()->getContainer("TpcDigiPar");
  TpcDigiMapper::getInstance()->init(tpcpar);
  //TpcDigiMapper::getInstance()->forceManualDriftVel(kTRUE);
  TpcAlignmentManager::init(tpcpar->getAlignmentFile());
  //TpcAlignmentManager::getInstance();
  cout<<"********************************************"<<endl;
  std::cout<<"pads:"<<TpcDigiMapper::getInstance()->getPadPlane()->GetNPads()<<std::endl;
  rtdb->print();
  std::cout<<"start:"<<nEvStart<<" end:"<<nEvents<<std::endl;
  //rtdb->saveOutput();

  fRun->Run(nEvStart,nEvents);

  timer.Stop();
  Double_t rtime = timer.RealTime();
  Double_t ctime = timer.CpuTime();
  printf("RealTime=%f seconds, CpuTime=%f seconds\n",rtime,ctime);
  
  std::cout<<"OutputFile: "<<outFile<<std::endl;
  
}