[u/mrichter/AliRoot.git] / TPC / Upgrade / macros / SCdistortionsForTDRfinal.C

void SCdistortionsForTDRfinalAll(){

  SCdistortionsForTDRfinalsingle(1,1,0);
  SCdistortionsForTDRfinalsingle(1,1,1);
  SCdistortionsForTDRfinalsingle(1,1,2);
  SCdistortionsForTDRfinalsingle(1,1,3);
  SCdistortionsForTDRfinalsingle(1,1,4);

}

void SCdistortionsForTDRfinal(Double_t radiusScale=1.5, Int_t epsScale=1,Int_t iOmegaTau = 1){
  //
  // do for one file (given by directory) the correction (specify the gas = iOmegaTau)
  // use the integrate along drift line option
  //
  // 1. Initialzation form space charge maps
  //
  AliTPCSpaceCharge3D *spaceCharge = new AliTPCSpaceCharge3D;
  const Double_t fgke0 = 8.854187817e-12; // vacuum permittivity [A·s/(V·m)]

  // omega tau parameters and TF1
  const Int_t nEps = 2;
  Int_t eps[nEps] = {20,10};
  Int_t col[nEps] = {kBlack,kBlack};

  const Int_t nOmegaTau = 5;
  Double_t omegaTau[nOmegaTau] = {0.34,0.32,0.43,1.77,1.84};
  TString tGas[nOmegaTau] = {"NeCO2","NeCO2_2","ArCO2","NeCF4","NeCF4_2"}; // CF4 is the same as CO2 here, but different omegaTau
  TString sGas[nOmegaTau] = {"Ne-CO_{2} (90-10)","Ne-CO_{2}-N_{2} (90-10-5)","Ar-CO_{2} (90-10)","Ne-CF_{4} (90-10)","Ne-CF_{4} (80-20)"};
  TF2 * fdiffR[nEps];
  TF2 * fdiffPhiR[nEps];
  TH2F * hdiffR[nEps];
  TH2F * hdiffPhiR[nEps];
  TF2 * fdiffIntR[nEps];
  TF2 * fdiffIntPhiR[nEps];
  TF2 * fdiffIntZ[nEps];
  TH2F * hdiffIntR[nEps];
  TH2F * hdiffIntPhiR[nEps];
  TH2F * hdiffIntZ[nEps];
  TH2F * hMap[nEps];
  TH2F * hDistRMap[nEps];
  TH2F * hDistRPMap[nEps];

  //use always the integrate option here
  Double_t integrateStep = 1.;

  TCanvas *cMap = new TCanvas("cMap","cMap",1200,500);
  cMap->Divide(2,1);

  TCanvas *cDistRMap = new TCanvas("cDistRMap","cDistRMap",1200,500);
  cDistRMap->Divide(2,1);

  TCanvas *cDistRPMap = new TCanvas("cDistRPMap","cDistRPMap",1200,500);
  cDistRPMap->Divide(2,1);

  TCanvas *cDistRNonIntMap = new TCanvas("cDistRNonIntMap","cDistRNonIntMap",1200,500);
  cDistRNonIntMap->Divide(2,1);

  TCanvas *cDistRPNonIntMap = new TCanvas("cDistRPNonIntMap","cDistRPNonIntMap",1200,500);
  cDistRPNonIntMap->Divide(2,1);

  TCanvas *cDistRIntMap = new TCanvas("cDistRIntMap","cDistRIntMap",1200,500);
  cDistRIntMap->Divide(2,1);

  TCanvas *cDistRPIntMap = new TCanvas("cDistRPIntMap","cDistRPIntMap",1200,500);
  cDistRPIntMap->Divide(2,1);

  TCanvas *cDistZIntMap = new TCanvas("cDistZIntMap","cDistZIntMap",1200,500);
  cDistZIntMap->Divide(2,1);


  TString outfilename = Form("SCdistortions_%s_50kHz_radiusScaling%.0f_epsScaling%d",tGas[iOmegaTau].Data(),radiusScale,epsScale);
  if(radiusScale>1.1 && radiusScale < 1.9) outfilename = Form("SCdistortions_%s_50kHz_radiusScaling%.1f_epsScaling%d",tGas[iOmegaTau].Data(),radiusScale,epsScale);

  //loop over epsilons
  for(Int_t iEps = 0; iEps < nEps; ++iEps){

    cMap->cd(iEps+1);

  // select gas 
  // 0 = NeCO2 
  // 1 = NeCO2N2 
  // 2 = ArCO2 
  // 3 = NeCF4 
  // 4 = ArCF4 
     if(radiusScale>1.1 && radiusScale < 1.9){
      cout<<"Open file = "<<Form("/Users/physics/ALICE/TPCupgrade/SpaceCharge/Maps/SC_%s_eps%d_50kHz_radiusScaling%.1f_epsScaling%d/SpaceChargeMap.root",tGas[iOmegaTau].Data(),eps[iEps],radiusScale,epsScale)<<endl;
      spaceCharge->SetSCDataFileName(Form("/Users/physics/ALICE/TPCupgrade/SpaceCharge/Maps/SC_%s_eps%d_50kHz_radiusScaling%.1f_epsScaling%d/SpaceChargeMap.root",tGas[iOmegaTau].Data(),eps[iEps],radiusScale,epsScale));
    }
    else{
      cout<<"Open file = "<<Form("/Users/physics/ALICE/TPCupgrade/SpaceCharge/Maps/SC_%s_eps%d_50kHz_radiusScaling%.0f_epsScaling%d/SpaceChargeMap.root",tGas[iOmegaTau].Data(),eps[iEps],radiusScale,epsScale)<<endl;
      spaceCharge->SetSCDataFileName(Form("/Users/physics/ALICE/TPCupgrade/SpaceCharge/Maps/SC_%s_eps%d_50kHz_radiusScaling%.0f_epsScaling%d/SpaceChargeMap.root",tGas[iOmegaTau].Data(),eps[iEps],radiusScale,epsScale));
    }

  // select omegaTau value
  if(iOmegaTau ==  1){
    spaceCharge->SetOmegaTauT1T2(omegaTau[iOmegaTau],1.00,0.99); // Ne CO2 N2 (90-10-5)
  }
  if(iOmegaTau ==  2){
    spaceCharge->SetOmegaTauT1T2(omegaTau[iOmegaTau],0.99,1.03); // Ar CO2 (90-10)
  }
  else if(iOmegaTau == 3){
    spaceCharge->SetOmegaTauT1T2(omegaTau[iOmegaTau],0.41,0.70); // Ne CF4 (90-10)
  }
  else if(iOmegaTau == 4){
    spaceCharge->SetOmegaTauT1T2(omegaTau[iOmegaTau],0.41,0.70); // Ne CF4 (80-20) (not in table use same as for other)
  }
  else{
    spaceCharge->SetOmegaTauT1T2(omegaTau[iOmegaTau],1,1.01); // Ne CO2 (90-10)
  }
  //
  //
  // init and add to corrections
  spaceCharge->InitSpaceCharge3DDistortion();
  spaceCharge->AddVisualCorrection(spaceCharge,1);

  // draw the map
  hMap[iEps] = (TH2F*)spaceCharge->CreateHistoSCinZR(0.);
  hMap[iEps]->Scale(fgke0/1e6*1e15); // C/m^3/e0 --> fC/cm^3
  hMap[iEps]->GetXaxis()->SetTitle("z (cm)");
  hMap[iEps]->GetYaxis()->SetTitle("r (cm)");
  hMap[iEps]->GetZaxis()->SetTitle("");
  //hMap[iEps]->GetZaxis()->SetTitle("#rho_{SC} (fC/cm^{3})");
  hMap[iEps]->SetTitleSize(0.05,"XYZ");
  //hMap[iEps]->SetTitleOffset(1.5,"XY");
  //hMap[iEps]->SetTitleOffset(0.9,"Z");
  hMap[iEps]->SetTitle(Form("#rho_{SC} (fC/cm^{3}) for %s, 50 kHz, #varepsilon = %d",sGas[iOmegaTau].Data(),eps[iEps]));
  hMap[iEps]->DrawCopy("colz");

  // draw the distortion maps
  cDistRMap->cd(iEps+1);
  hDistRMap[iEps] = (TH2F*)spaceCharge->CreateHistoDRinZR(0.);
  hDistRMap[iEps]->SetMaximum(25);
  hDistRMap[iEps]->SetMinimum(-12);
  hDistRMap[iEps]->GetXaxis()->SetTitle("z (cm)");
  hDistRMap[iEps]->GetYaxis()->SetTitle("r (cm)");
  hDistRMap[iEps]->GetZaxis()->SetTitle("");
  hDistRMap[iEps]->SetTitleSize(0.05,"XYZ");
  //hDistRMap[iEps]->SetTitleOffset(1.5,"XY");
  //hDistRMap[iEps]->SetTitleOffset(0.9,"Z");
  hDistRMap[iEps]->SetTitle(Form("dr (cm) for %s, 50 kHz, #varepsilon = %d",sGas[iOmegaTau].Data(),eps[iEps]));
  hDistRMap[iEps]->DrawCopy("colz");

  cDistRPMap->cd(iEps+1);
  hDistRPMap[iEps] = (TH2F*)spaceCharge->CreateHistoDRinZR(0.);
  hDistRPMap[iEps]->SetMaximum(25);
  hDistRPMap[iEps]->SetMinimum(-12);
  hDistRPMap[iEps]->GetXaxis()->SetTitle("z (cm)");
  hDistRPMap[iEps]->GetYaxis()->SetTitle("r (cm)");
  hDistRPMap[iEps]->GetZaxis()->SetTitle("");
  hDistRPMap[iEps]->SetTitleSize(0.05,"XYZ");
  //hDistRPMap[iEps]->SetTitleOffset(1.5,"XY");
  //hDistRPMap[iEps]->SetTitleOffset(0.9,"Z");
  hDistRPMap[iEps]->SetTitle(Form("d(r#varphi) (cm) for %s, 50 kHz, #varepsilon = %d",sGas[iOmegaTau].Data(),eps[iEps]));
  hDistRPMap[iEps]->DrawCopy("colz");

  
  //
  // 2. get TF2 with differences 
  //
  // get corrections (at y = 0) for visual correction 1 (last argument)
  // 0 = dR
  // 1 = dPhiR
  // 2 = dZ?    
  fdiffR[iEps]       = new TF2(Form("fdiffR%d",iEps), Form("AliTPCCorrection::GetDistXYZ(y,0,x,0,1)"),-250,250,85,250);
  fdiffPhiR[iEps]    = new TF2(Form("fdiffPhiR%d",iEps), Form("AliTPCCorrection::GetDistXYZ(y,0,x,1,1)"),-250,250,85,250);
  hdiffR[iEps] = (TH2F*)fdiffR[iEps]->GetHistogram();
  hdiffPhiR[iEps] = (TH2F*)fdiffPhiR[iEps]->GetHistogram();

  hdiffR[iEps]->SetName(fdiffR[iEps]->GetName());
  hdiffPhiR[iEps]->SetName(fdiffPhiR[iEps]->GetName());

  fdiffIntR[iEps]       = new TF2(Form("fdiffIntR%d",iEps), Form("AliTPCCorrection::GetDistXYZIntegrateZ(y,0,x,0,1,%f)",integrateStep),-250,250,85,250);
  fdiffIntPhiR[iEps]    = new TF2(Form("fdiffIntPhiR%d",iEps), Form("AliTPCCorrection::GetDistXYZIntegrateZ(y,0,x,1,1,%f)",integrateStep),-250,250,85,250);
  fdiffIntZ[iEps]    = new TF2(Form("fdiffIntZ%d",iEps), Form("AliTPCCorrection::GetDistXYZIntegrateZ(y,0,x,2,1,%f)",integrateStep),-250,250,85,250);
  hdiffIntR[iEps] = (TH2F*)fdiffIntR[iEps]->GetHistogram();
  hdiffIntPhiR[iEps] = (TH2F*)fdiffIntPhiR[iEps]->GetHistogram();
  hdiffIntZ[iEps] = (TH2F*)fdiffIntZ[iEps]->GetHistogram();

  hdiffIntR[iEps]->SetName(fdiffIntR[iEps]->GetName());
  hdiffIntPhiR[iEps]->SetName(fdiffIntPhiR[iEps]->GetName());
  hdiffIntZ[iEps]->SetName(fdiffIntZ[iEps]->GetName());

  
  //
  // 3. Plot and store TH1Fs
  //
  
  // draw the distortion maps
  cDistRNonIntMap->cd(iEps+1);
  hdiffR[iEps]->SetMaximum(25);
  hdiffR[iEps]->SetMinimum(-12);
  hdiffR[iEps]->GetXaxis()->SetTitle("z (cm)");
  hdiffR[iEps]->GetYaxis()->SetTitle("r (cm)");
  //hdiffR[iEps]->GetZaxis()->SetTitle("dr (cm)");
  hdiffR[iEps]->SetTitleSize(0.05,"XYZ");
  //hdiffR[iEps]->SetTitleOffset(1.5,"XY");
  //hdiffR[iEps]->SetTitleOffset(0.9,"Z");
  hdiffR[iEps]->SetTitle(Form("dr (cm) for %s, 50 kHz, #varepsilon = %d",sGas[iOmegaTau].Data(),eps[iEps]));
  hdiffR[iEps]->DrawCopy("colz");


  cDistRPNonIntMap->cd(iEps+1);
  hdiffPhiR[iEps]->SetMaximum(25);
  hdiffPhiR[iEps]->SetMinimum(-12);
  hdiffPhiR[iEps]->GetXaxis()->SetTitle("z (cm)");
  hdiffPhiR[iEps]->GetYaxis()->SetTitle("r (cm)");
  hdiffPhiR[iEps]->SetTitleSize(0.05,"XYZ");
  //hdiffPhiR[iEps]->SetTitleOffset(1.5,"XY");
  //hdiffPhiR[iEps]->SetTitleOffset(0.9,"Z");
  hdiffPhiR[iEps]->SetTitle(Form("d(r#varphi) (cm) for %s, 50 kHz, #varepsilon = %d",sGas[iOmegaTau].Data(),eps[iEps]));
  hdiffPhiR[iEps]->DrawCopy("colz");

  cDistRIntMap->cd(iEps+1);
  hdiffIntR[iEps]->SetMaximum(25);
  hdiffIntR[iEps]->SetMinimum(-12);
  hdiffIntR[iEps]->GetXaxis()->SetTitle("z (cm)");
  hdiffIntR[iEps]->GetYaxis()->SetTitle("r (cm)");
  //hdiffIntR[iEps]->GetZaxis()->SetTitle("dr (cm)");
  hdiffIntR[iEps]->SetTitleSize(0.05,"XYZ");
  //hdiffIntR[iEps]->SetTitleOffset(1.5,"XY");
  //hdiffIntR[iEps]->SetTitleOffset(0.9,"Z");
  hdiffIntR[iEps]->SetTitle(Form("dr (cm) for %s, 50 kHz, #varepsilon = %d",sGas[iOmegaTau].Data(),eps[iEps]));
  hdiffIntR[iEps]->DrawCopy("colz");


  cDistRPIntMap->cd(iEps+1);
  hdiffIntPhiR[iEps]->SetMaximum(5);
  hdiffIntPhiR[iEps]->SetMinimum(-9);
  hdiffIntPhiR[iEps]->GetXaxis()->SetTitle("z (cm)");
  hdiffIntPhiR[iEps]->GetYaxis()->SetTitle("r (cm)");
  hdiffIntPhiR[iEps]->SetTitleSize(0.05,"XYZ");
  //hdiffIntPhiR[iEps]->SetTitleOffset(1.5,"XY");
  //hdiffIntPhiR[iEps]->SetTitleOffset(0.9,"Z");
  hdiffIntPhiR[iEps]->SetTitle(Form("d(r#varphi) (cm) for %s, 50 kHz, #varepsilon = %d",sGas[iOmegaTau].Data(),eps[iEps]));
  hdiffIntPhiR[iEps]->DrawCopy("colz");

  cDistZIntMap->cd(iEps+1);
  hdiffIntZ[iEps]->SetMaximum(5);
  hdiffIntZ[iEps]->SetMinimum(-5);
  hdiffIntZ[iEps]->GetXaxis()->SetTitle("z (cm)");
  hdiffIntZ[iEps]->GetYaxis()->SetTitle("r (cm)");
  hdiffIntZ[iEps]->SetTitleSize(0.05,"XYZ");
  //hdiffIntZ[iEps]->SetTitleOffset(1.5,"XY");
  //hdiffIntZ[iEps]->SetTitleOffset(0.9,"Z");
  hdiffIntZ[iEps]->SetTitle(Form("dz (cm) for %s, 50 kHz, #varepsilon = %d",sGas[iOmegaTau].Data(),eps[iEps]));
  hdiffIntZ[iEps]->DrawCopy("colz");
  }
  
  // just control histograms (not saved)
  // cDistRNonIntMap->SaveAs(Form("%s_TDR_DistortR.eps",outfilename.Data()));
  // cDistRNonIntMap->SaveAs(Form("%s_TDR_DistortR.png",outfilename.Data()));
  // cDistRNonIntMap->SaveAs(Form("%s_TDR_DistortR.pdf",outfilename.Data()));
  // cDistRPNonIntMap->SaveAs(Form("%s_TDR_DistortRPhi.eps",outfilename.Data()));
  // cDistRPNonIntMap->SaveAs(Form("%s_TDR_DistortRPhi.png",outfilename.Data()));
  // cDistRPNonIntMap->SaveAs(Form("%s_TDR_DistortRPhi.pdf",outfilename.Data()));

  // histograms for TDR (saved)
  cMap->SaveAs(Form("%s_TDR_SpaceCharge.eps",outfilename.Data()));
  cMap->SaveAs(Form("%s_TDR_SpaceCharge.png",outfilename.Data()));
  cMap->SaveAs(Form("%s_TDR_SpaceCharge.pdf",outfilename.Data()));
  cDistRIntMap->SaveAs(Form("%s_TDR_DistortIntR.eps",outfilename.Data()));
  cDistRIntMap->SaveAs(Form("%s_TDR_DistortIntR.png",outfilename.Data()));
  cDistRIntMap->SaveAs(Form("%s_TDR_DistortIntR.pdf",outfilename.Data()));
  cDistRPIntMap->SaveAs(Form("%s_TDR_DistortIntRPhi.eps",outfilename.Data()));
  cDistRPIntMap->SaveAs(Form("%s_TDR_DistortIntRPhi.png",outfilename.Data()));
  cDistRPIntMap->SaveAs(Form("%s_TDR_DistortIntRPhi.pdf",outfilename.Data()));
  cDistZIntMap->SaveAs(Form("%s_TDR_DistortIntZ.eps",outfilename.Data()));
  cDistZIntMap->SaveAs(Form("%s_TDR_DistortIntZ.png",outfilename.Data()));
  cDistZIntMap->SaveAs(Form("%s_TDR_DistortIntZ.pdf",outfilename.Data()));
}

  
//____________________________________________________________//
void setupLegend(TLegend *currentLegend=0,float currentTextSize=0.07){
  currentLegend->SetTextFont(42);
  currentLegend->SetBorderSize(0);
  currentLegend->SetFillStyle(0);
  currentLegend->SetFillColor(0);
  currentLegend->SetMargin(0.25);
  currentLegend->SetTextSize(currentTextSize);
  currentLegend->SetEntrySeparation(0.5);
  return;
}
Commit	Line	Data
9cfdb8b6	1	void SCdistortionsForTDRfinalAll(){
	2
	3	SCdistortionsForTDRfinalsingle(1,1,0);
	4	SCdistortionsForTDRfinalsingle(1,1,1);
	5	SCdistortionsForTDRfinalsingle(1,1,2);
	6	SCdistortionsForTDRfinalsingle(1,1,3);
	7	SCdistortionsForTDRfinalsingle(1,1,4);
	8
	9	}
	10
	11	void SCdistortionsForTDRfinal(Double_t radiusScale=1.5, Int_t epsScale=1,Int_t iOmegaTau = 1){
	12	//
	13	// do for one file (given by directory) the correction (specify the gas = iOmegaTau)
	14	// use the integrate along drift line option
	15	//
	16	// 1. Initialzation form space charge maps
	17	//
	18	AliTPCSpaceCharge3D *spaceCharge = new AliTPCSpaceCharge3D;
	19	const Double_t fgke0 = 8.854187817e-12; // vacuum permittivity [A·s/(V·m)]
	20
	21	// omega tau parameters and TF1
	22	const Int_t nEps = 2;
	23	Int_t eps[nEps] = {20,10};
	24	Int_t col[nEps] = {kBlack,kBlack};
	25
	26	const Int_t nOmegaTau = 5;
	27	Double_t omegaTau[nOmegaTau] = {0.34,0.32,0.43,1.77,1.84};
	28	TString tGas[nOmegaTau] = {"NeCO2","NeCO2_2","ArCO2","NeCF4","NeCF4_2"}; // CF4 is the same as CO2 here, but different omegaTau
21394c38	29	TString sGas[nOmegaTau] = {"Ne-CO_{2} (90-10)","Ne-CO_{2}-N_{2} (90-10-5)","Ar-CO_{2} (90-10)","Ne-CF_{4} (90-10)","Ne-CF_{4} (80-20)"};
9cfdb8b6	30	TF2 * fdiffR[nEps];
	31	TF2 * fdiffPhiR[nEps];
	32	TH2F * hdiffR[nEps];
	33	TH2F * hdiffPhiR[nEps];
	34	TF2 * fdiffIntR[nEps];
	35	TF2 * fdiffIntPhiR[nEps];
	36	TF2 * fdiffIntZ[nEps];
	37	TH2F * hdiffIntR[nEps];
	38	TH2F * hdiffIntPhiR[nEps];
	39	TH2F * hdiffIntZ[nEps];
	40	TH2F * hMap[nEps];
	41	TH2F * hDistRMap[nEps];
	42	TH2F * hDistRPMap[nEps];
	43
	44	//use always the integrate option here
	45	Double_t integrateStep = 1.;
	46
	47	TCanvas *cMap = new TCanvas("cMap","cMap",1200,500);
	48	cMap->Divide(2,1);
	49
	50	TCanvas *cDistRMap = new TCanvas("cDistRMap","cDistRMap",1200,500);
	51	cDistRMap->Divide(2,1);
	52
	53	TCanvas *cDistRPMap = new TCanvas("cDistRPMap","cDistRPMap",1200,500);
	54	cDistRPMap->Divide(2,1);
	55
	56	TCanvas *cDistRNonIntMap = new TCanvas("cDistRNonIntMap","cDistRNonIntMap",1200,500);
	57	cDistRNonIntMap->Divide(2,1);
	58
	59	TCanvas *cDistRPNonIntMap = new TCanvas("cDistRPNonIntMap","cDistRPNonIntMap",1200,500);
	60	cDistRPNonIntMap->Divide(2,1);
	61
	62	TCanvas *cDistRIntMap = new TCanvas("cDistRIntMap","cDistRIntMap",1200,500);
	63	cDistRIntMap->Divide(2,1);
	64
	65	TCanvas *cDistRPIntMap = new TCanvas("cDistRPIntMap","cDistRPIntMap",1200,500);
	66	cDistRPIntMap->Divide(2,1);
	67
	68	TCanvas *cDistZIntMap = new TCanvas("cDistZIntMap","cDistZIntMap",1200,500);
	69	cDistZIntMap->Divide(2,1);
	70
	71
	72
	73	TString outfilename = Form("SCdistortions_%s_50kHz_radiusScaling%.0f_epsScaling%d",tGas[iOmegaTau].Data(),radiusScale,epsScale);
	74	if(radiusScale>1.1 && radiusScale < 1.9) outfilename = Form("SCdistortions_%s_50kHz_radiusScaling%.1f_epsScaling%d",tGas[iOmegaTau].Data(),radiusScale,epsScale);
	75
	76	//loop over epsilons
	77	for(Int_t iEps = 0; iEps < nEps; ++iEps){
	78
	79	cMap->cd(iEps+1);
	80
	81	// select gas
	82	// 0 = NeCO2
	83	// 1 = NeCO2N2
	84	// 2 = ArCO2
	85	// 3 = NeCF4
	86	// 4 = ArCF4
	87	if(radiusScale>1.1 && radiusScale < 1.9){
	88	cout<<"Open file = "<<Form("/Users/physics/ALICE/TPCupgrade/SpaceCharge/Maps/SC_%s_eps%d_50kHz_radiusScaling%.1f_epsScaling%d/SpaceChargeMap.root",tGas[iOmegaTau].Data(),eps[iEps],radiusScale,epsScale)<<endl;
	89	spaceCharge->SetSCDataFileName(Form("/Users/physics/ALICE/TPCupgrade/SpaceCharge/Maps/SC_%s_eps%d_50kHz_radiusScaling%.1f_epsScaling%d/SpaceChargeMap.root",tGas[iOmegaTau].Data(),eps[iEps],radiusScale,epsScale));
	90	}
	91	else{
	92	cout<<"Open file = "<<Form("/Users/physics/ALICE/TPCupgrade/SpaceCharge/Maps/SC_%s_eps%d_50kHz_radiusScaling%.0f_epsScaling%d/SpaceChargeMap.root",tGas[iOmegaTau].Data(),eps[iEps],radiusScale,epsScale)<<endl;
	93	spaceCharge->SetSCDataFileName(Form("/Users/physics/ALICE/TPCupgrade/SpaceCharge/Maps/SC_%s_eps%d_50kHz_radiusScaling%.0f_epsScaling%d/SpaceChargeMap.root",tGas[iOmegaTau].Data(),eps[iEps],radiusScale,epsScale));
94	}
95
96	// select omegaTau value
97	if(iOmegaTau == 1){
98	spaceCharge->SetOmegaTauT1T2(omegaTau[iOmegaTau],1.00,0.99); // Ne CO2 N2 (90-10-5)
99	}
100	if(iOmegaTau == 2){
101	spaceCharge->SetOmegaTauT1T2(omegaTau[iOmegaTau],0.99,1.03); // Ar CO2 (90-10)
102	}
103	else if(iOmegaTau == 3){
104	spaceCharge->SetOmegaTauT1T2(omegaTau[iOmegaTau],0.41,0.70); // Ne CF4 (90-10)
105	}
106	else if(iOmegaTau == 4){
107	spaceCharge->SetOmegaTauT1T2(omegaTau[iOmegaTau],0.41,0.70); // Ne CF4 (80-20) (not in table use same as for other)
108	}
109	else{
110	spaceCharge->SetOmegaTauT1T2(omegaTau[iOmegaTau],1,1.01); // Ne CO2 (90-10)
111	}
112	//
113	//
114	// init and add to corrections
115	spaceCharge->InitSpaceCharge3DDistortion();
116	spaceCharge->AddVisualCorrection(spaceCharge,1);
117
118	// draw the map
119	hMap[iEps] = (TH2F*)spaceCharge->CreateHistoSCinZR(0.);
120	hMap[iEps]->Scale(fgke0/1e6*1e15); // C/m^3/e0 --> fC/cm^3
121	hMap[iEps]->GetXaxis()->SetTitle("z (cm)");
122	hMap[iEps]->GetYaxis()->SetTitle("r (cm)");
123	hMap[iEps]->GetZaxis()->SetTitle("");
124	//hMap[iEps]->GetZaxis()->SetTitle("#rho_{SC} (fC/cm^{3})");
125	hMap[iEps]->SetTitleSize(0.05,"XYZ");
126	//hMap[iEps]->SetTitleOffset(1.5,"XY");
127	//hMap[iEps]->SetTitleOffset(0.9,"Z");
21394c38	128	hMap[iEps]->SetTitle(Form("#rho_{SC} (fC/cm^{3}) for %s, 50 kHz, #varepsilon = %d",sGas[iOmegaTau].Data(),eps[iEps]));
9cfdb8b6	129	hMap[iEps]->DrawCopy("colz");
	130
	131	// draw the distortion maps
	132	cDistRMap->cd(iEps+1);
	133	hDistRMap[iEps] = (TH2F*)spaceCharge->CreateHistoDRinZR(0.);
	134	hDistRMap[iEps]->SetMaximum(25);
	135	hDistRMap[iEps]->SetMinimum(-12);
	136	hDistRMap[iEps]->GetXaxis()->SetTitle("z (cm)");
	137	hDistRMap[iEps]->GetYaxis()->SetTitle("r (cm)");
	138	hDistRMap[iEps]->GetZaxis()->SetTitle("");
	139	hDistRMap[iEps]->SetTitleSize(0.05,"XYZ");
	140	//hDistRMap[iEps]->SetTitleOffset(1.5,"XY");
	141	//hDistRMap[iEps]->SetTitleOffset(0.9,"Z");
21394c38	142	hDistRMap[iEps]->SetTitle(Form("dr (cm) for %s, 50 kHz, #varepsilon = %d",sGas[iOmegaTau].Data(),eps[iEps]));
9cfdb8b6	143	hDistRMap[iEps]->DrawCopy("colz");
	144
	145	cDistRPMap->cd(iEps+1);
	146	hDistRPMap[iEps] = (TH2F*)spaceCharge->CreateHistoDRinZR(0.);
	147	hDistRPMap[iEps]->SetMaximum(25);
	148	hDistRPMap[iEps]->SetMinimum(-12);
	149	hDistRPMap[iEps]->GetXaxis()->SetTitle("z (cm)");
	150	hDistRPMap[iEps]->GetYaxis()->SetTitle("r (cm)");
	151	hDistRPMap[iEps]->GetZaxis()->SetTitle("");
	152	hDistRPMap[iEps]->SetTitleSize(0.05,"XYZ");
	153	//hDistRPMap[iEps]->SetTitleOffset(1.5,"XY");
	154	//hDistRPMap[iEps]->SetTitleOffset(0.9,"Z");
21394c38	155	hDistRPMap[iEps]->SetTitle(Form("d(r#varphi) (cm) for %s, 50 kHz, #varepsilon = %d",sGas[iOmegaTau].Data(),eps[iEps]));
9cfdb8b6	156	hDistRPMap[iEps]->DrawCopy("colz");
	157
	158
	159	//
	160	// 2. get TF2 with differences
	161	//
	162	// get corrections (at y = 0) for visual correction 1 (last argument)
	163	// 0 = dR
	164	// 1 = dPhiR
	165	// 2 = dZ?
	166	fdiffR[iEps] = new TF2(Form("fdiffR%d",iEps), Form("AliTPCCorrection::GetDistXYZ(y,0,x,0,1)"),-250,250,85,250);
	167	fdiffPhiR[iEps] = new TF2(Form("fdiffPhiR%d",iEps), Form("AliTPCCorrection::GetDistXYZ(y,0,x,1,1)"),-250,250,85,250);
	168	hdiffR[iEps] = (TH2F*)fdiffR[iEps]->GetHistogram();
	169	hdiffPhiR[iEps] = (TH2F*)fdiffPhiR[iEps]->GetHistogram();
	170
	171	hdiffR[iEps]->SetName(fdiffR[iEps]->GetName());
	172	hdiffPhiR[iEps]->SetName(fdiffPhiR[iEps]->GetName());
	173
	174	fdiffIntR[iEps] = new TF2(Form("fdiffIntR%d",iEps), Form("AliTPCCorrection::GetDistXYZIntegrateZ(y,0,x,0,1,%f)",integrateStep),-250,250,85,250);
	175	fdiffIntPhiR[iEps] = new TF2(Form("fdiffIntPhiR%d",iEps), Form("AliTPCCorrection::GetDistXYZIntegrateZ(y,0,x,1,1,%f)",integrateStep),-250,250,85,250);
	176	fdiffIntZ[iEps] = new TF2(Form("fdiffIntZ%d",iEps), Form("AliTPCCorrection::GetDistXYZIntegrateZ(y,0,x,2,1,%f)",integrateStep),-250,250,85,250);
	177	hdiffIntR[iEps] = (TH2F*)fdiffIntR[iEps]->GetHistogram();
	178	hdiffIntPhiR[iEps] = (TH2F*)fdiffIntPhiR[iEps]->GetHistogram();
	179	hdiffIntZ[iEps] = (TH2F*)fdiffIntZ[iEps]->GetHistogram();
	180
	181	hdiffIntR[iEps]->SetName(fdiffIntR[iEps]->GetName());
	182	hdiffIntPhiR[iEps]->SetName(fdiffIntPhiR[iEps]->GetName());
	183	hdiffIntZ[iEps]->SetName(fdiffIntZ[iEps]->GetName());
	184
	185
	186
	187	//
	188	// 3. Plot and store TH1Fs
	189	//
	190
	191	// draw the distortion maps
	192	cDistRNonIntMap->cd(iEps+1);
	193	hdiffR[iEps]->SetMaximum(25);
	194	hdiffR[iEps]->SetMinimum(-12);
	195	hdiffR[iEps]->GetXaxis()->SetTitle("z (cm)");
	196	hdiffR[iEps]->GetYaxis()->SetTitle("r (cm)");
	197	//hdiffR[iEps]->GetZaxis()->SetTitle("dr (cm)");
	198	hdiffR[iEps]->SetTitleSize(0.05,"XYZ");
	199	//hdiffR[iEps]->SetTitleOffset(1.5,"XY");
	200	//hdiffR[iEps]->SetTitleOffset(0.9,"Z");
21394c38	201	hdiffR[iEps]->SetTitle(Form("dr (cm) for %s, 50 kHz, #varepsilon = %d",sGas[iOmegaTau].Data(),eps[iEps]));
9cfdb8b6	202	hdiffR[iEps]->DrawCopy("colz");
	203
	204
	205	cDistRPNonIntMap->cd(iEps+1);
	206	hdiffPhiR[iEps]->SetMaximum(25);
	207	hdiffPhiR[iEps]->SetMinimum(-12);
	208	hdiffPhiR[iEps]->GetXaxis()->SetTitle("z (cm)");
	209	hdiffPhiR[iEps]->GetYaxis()->SetTitle("r (cm)");
	210	hdiffPhiR[iEps]->SetTitleSize(0.05,"XYZ");
	211	//hdiffPhiR[iEps]->SetTitleOffset(1.5,"XY");
	212	//hdiffPhiR[iEps]->SetTitleOffset(0.9,"Z");
21394c38	213	hdiffPhiR[iEps]->SetTitle(Form("d(r#varphi) (cm) for %s, 50 kHz, #varepsilon = %d",sGas[iOmegaTau].Data(),eps[iEps]));
9cfdb8b6	214	hdiffPhiR[iEps]->DrawCopy("colz");
	215
	216	cDistRIntMap->cd(iEps+1);
	217	hdiffIntR[iEps]->SetMaximum(25);
	218	hdiffIntR[iEps]->SetMinimum(-12);
	219	hdiffIntR[iEps]->GetXaxis()->SetTitle("z (cm)");
	220	hdiffIntR[iEps]->GetYaxis()->SetTitle("r (cm)");
	221	//hdiffIntR[iEps]->GetZaxis()->SetTitle("dr (cm)");
	222	hdiffIntR[iEps]->SetTitleSize(0.05,"XYZ");
	223	//hdiffIntR[iEps]->SetTitleOffset(1.5,"XY");
	224	//hdiffIntR[iEps]->SetTitleOffset(0.9,"Z");
21394c38	225	hdiffIntR[iEps]->SetTitle(Form("dr (cm) for %s, 50 kHz, #varepsilon = %d",sGas[iOmegaTau].Data(),eps[iEps]));
9cfdb8b6	226	hdiffIntR[iEps]->DrawCopy("colz");
	227
	228
	229	cDistRPIntMap->cd(iEps+1);
	230	hdiffIntPhiR[iEps]->SetMaximum(5);
	231	hdiffIntPhiR[iEps]->SetMinimum(-9);
	232	hdiffIntPhiR[iEps]->GetXaxis()->SetTitle("z (cm)");
	233	hdiffIntPhiR[iEps]->GetYaxis()->SetTitle("r (cm)");
	234	hdiffIntPhiR[iEps]->SetTitleSize(0.05,"XYZ");
	235	//hdiffIntPhiR[iEps]->SetTitleOffset(1.5,"XY");
	236	//hdiffIntPhiR[iEps]->SetTitleOffset(0.9,"Z");
21394c38	237	hdiffIntPhiR[iEps]->SetTitle(Form("d(r#varphi) (cm) for %s, 50 kHz, #varepsilon = %d",sGas[iOmegaTau].Data(),eps[iEps]));
9cfdb8b6	238	hdiffIntPhiR[iEps]->DrawCopy("colz");
	239
	240	cDistZIntMap->cd(iEps+1);
	241	hdiffIntZ[iEps]->SetMaximum(5);
	242	hdiffIntZ[iEps]->SetMinimum(-5);
	243	hdiffIntZ[iEps]->GetXaxis()->SetTitle("z (cm)");
	244	hdiffIntZ[iEps]->GetYaxis()->SetTitle("r (cm)");
	245	hdiffIntZ[iEps]->SetTitleSize(0.05,"XYZ");
	246	//hdiffIntZ[iEps]->SetTitleOffset(1.5,"XY");
	247	//hdiffIntZ[iEps]->SetTitleOffset(0.9,"Z");
21394c38	248	hdiffIntZ[iEps]->SetTitle(Form("dz (cm) for %s, 50 kHz, #varepsilon = %d",sGas[iOmegaTau].Data(),eps[iEps]));
9cfdb8b6	249	hdiffIntZ[iEps]->DrawCopy("colz");
	250	}
	251
	252	// just control histograms (not saved)
	253	// cDistRNonIntMap->SaveAs(Form("%s_TDR_DistortR.eps",outfilename.Data()));
	254	// cDistRNonIntMap->SaveAs(Form("%s_TDR_DistortR.png",outfilename.Data()));
	255	// cDistRNonIntMap->SaveAs(Form("%s_TDR_DistortR.pdf",outfilename.Data()));
	256	// cDistRPNonIntMap->SaveAs(Form("%s_TDR_DistortRPhi.eps",outfilename.Data()));
	257	// cDistRPNonIntMap->SaveAs(Form("%s_TDR_DistortRPhi.png",outfilename.Data()));
	258	// cDistRPNonIntMap->SaveAs(Form("%s_TDR_DistortRPhi.pdf",outfilename.Data()));
	259
	260	// histograms for TDR (saved)
	261	cMap->SaveAs(Form("%s_TDR_SpaceCharge.eps",outfilename.Data()));
	262	cMap->SaveAs(Form("%s_TDR_SpaceCharge.png",outfilename.Data()));
	263	cMap->SaveAs(Form("%s_TDR_SpaceCharge.pdf",outfilename.Data()));
	264	cDistRIntMap->SaveAs(Form("%s_TDR_DistortIntR.eps",outfilename.Data()));
	265	cDistRIntMap->SaveAs(Form("%s_TDR_DistortIntR.png",outfilename.Data()));
	266	cDistRIntMap->SaveAs(Form("%s_TDR_DistortIntR.pdf",outfilename.Data()));
	267	cDistRPIntMap->SaveAs(Form("%s_TDR_DistortIntRPhi.eps",outfilename.Data()));
	268	cDistRPIntMap->SaveAs(Form("%s_TDR_DistortIntRPhi.png",outfilename.Data()));
	269	cDistRPIntMap->SaveAs(Form("%s_TDR_DistortIntRPhi.pdf",outfilename.Data()));
	270	cDistZIntMap->SaveAs(Form("%s_TDR_DistortIntZ.eps",outfilename.Data()));
	271	cDistZIntMap->SaveAs(Form("%s_TDR_DistortIntZ.png",outfilename.Data()));
	272	cDistZIntMap->SaveAs(Form("%s_TDR_DistortIntZ.pdf",outfilename.Data()));
	273	}
	274
	275
	276	//____________________________________________________________//
	277	void setupLegend(TLegend *currentLegend=0,float currentTextSize=0.07){
	278	currentLegend->SetTextFont(42);
	279	currentLegend->SetBorderSize(0);
	280	currentLegend->SetFillStyle(0);
	281	currentLegend->SetFillColor(0);
	282	currentLegend->SetMargin(0.25);
	283	currentLegend->SetTextSize(currentTextSize);
	284	currentLegend->SetEntrySeparation(0.5);
	285	return;
	286	}