//ROOT script to create the deviation Map using Runge Kutta 
//Argument rBin forces TpcInhFieldDrifter to calculate only one row in r
//In this way the macro can effectively be run on the cluster
//Output will be rBin ascii files which need to be combined to one 
//generic Deviation Map file.

//Author: Felix Boehmer


#include <fstream>
#include <vector>
#include <string>
#include <iostream>
#include "TpcInhFieldDrifter.h"
#include "PndMultiField.h"
#include "PndTransMap.h"
#include "PndDipoleMap.h"
#include "PndSolenoidMap.h"
#include "TpcEFieldCyl.h"
#include "FairRunSim.h"


void createDevMapMulti(TString efieldfile,
		       TString paramFile,
		       int rBin,
		       TString outFile = "") {

  //define variables
  int rBins = 26;   
  int zBins = 149;
  int _rBin = rBin;

  FairRunSim *fRun = new FairRunSim();
  fRun->SetBeamMom(3.672);
  

  std::string basedir = "/afs/e18/panda/SIM/fboehmer/trunk/";
  //std::string targetdir(targetDir.Data());
  
  std::string efield_file(efieldfile.Data());
  
  // Use TpcEFieldCyl for building the B-Field (custom input) ------------
  
  std::string bfield_file = "andrea_b_field.dat"; //for bfield as efield
  // ------------------------------------------------------------------------
  TString out = outFile;
  std::string test(out.Data());
  //if no outfile was given, generate name out of efieldfile and rBin:
  if (test.length()<2) {
    out = efieldfile;
    char binBuff[5];
    sprintf(binBuff, "%i", rBin);
    TString fill("DevMap_rBin");
    fill.Append(binBuff);
    out.ReplaceAll("Efield", fill);
  }
  std::string out_file(out.Data());
  //std::string param_file("tpc/parfiles/tpc_edged.par");
  std::string param_file(paramFile.Data());
  
  
  // PndMultiField part -----------------------------------------------------
  
  PndMultiField *bField= new PndMultiField("FULL");
  
  // PndTransMap *map_tr= new PndTransMap("TransMap", "R");
  // PndDipoleMap *map_d1= new PndDipoleMap("DipoleMap", "R");
  // //  PndDipoleMap *map_d2= new PndDipoleMap("DipoleMap2", "R");
  // PndSolenoidMap *map_s1= new PndSolenoidMap("SolenoidMap1", "R");
  // PndSolenoidMap *map_s2= new PndSolenoidMap("SolenoidMap2", "R");
  // PndSolenoidMap *map_s3= new PndSolenoidMap("SolenoidMap3", "R");
  // PndSolenoidMap *map_s4= new PndSolenoidMap("SolenoidMap4", "R");
  
  // bField->AddField(map_tr);
  // bField->AddField(map_d1);
  // //bField->AddField(map_d2);
  // bField->AddField(map_s1);
  // bField->AddField(map_s2);
  // bField->AddField(map_s3);
  // bField->AddField(map_s4);
  
  bField->Init();
  
  //create Efield
  TpcEFieldCyl* eField = new TpcEFieldCyl(efield_file.c_str());
   
  //if(bField != 0)
  //  std::cout<<"\nsuccessfully initialized the B-Field";
  
  // ------------------------------------------------------------------------
  
  
  
  // start the simulation
  // TpcInhFieldDrifter* drifter = new TpcInhFieldDrifter(efield_file.c_str(),
  // 							   bfield_file.c_str(),
// 							   out_file.c_str(),
// 							   param_file.c_str(),
// 							   rBins, zBins, _rBin);
  
  TpcInhFieldDrifter* drifter = new TpcInhFieldDrifter(eField,
						       bField,
						       out_file.c_str(),
						       param_file.c_str(),
						       rBins, zBins, _rBin);
  
  
  std::cout<<"\nStarting the Drifter..."<<std::endl;
  
  drifter->run();
  delete drifter;
  //delete bField;

}