[u/mrichter/AliRoot.git] / TPC / CalibMacros / AliTPCCorrectionDemo.C

/*
 gSystem->SetIncludePath("-I$ROOTSYS/include -I$ALICE_ROOT/ -I$ALICE_ROOT/include -I$ALICE_ROOT/STEER -I$ALICE_ROOT/TPC -I$ALICE_ROOT/ITS -I$ALICE_ROOT/TRD -I$ALICE_ROOT/TOF -I$ALICE_ROOT/RAW -I$ALICE_ROOT/STAT");

 */
#if !defined(__CINT__) || defined(__MAKECINT__)
#include "THnSparse.h"
#include "TLatex.h"
#include "TCanvas.h"
#include "TLegend.h"
#include "TSystem.h"
#include "TFile.h"
#include "TChain.h"
#include "TCut.h"
#include "TH3.h"
#include "TH2F.h"
#include "TProfile3D.h"
#include "TMath.h" 
#include "TVectorD.h"
#include "TMatrixD.h"
#include "TStatToolkit.h"
#include "TTreeStream.h"
#include "AliExternalTrackParam.h"
#include "AliESDfriend.h"
#include "AliTPCcalibTime.h"
#include "TROOT.h"
#include "AliXRDPROOFtoolkit.h"
#include "AliTPCCorrection.h"
#include "AliTPCExBTwist.h"
#include "AliTPCGGVoltError.h"
#include "AliTPCComposedCorrection.h"
#include "AliTPCExBConical.h"
#include "TPostScript.h"
#include "TStyle.h"
#include "AliTrackerBase.h"
#include "TGraph.h"
#include "AliCDBManager.h"
#include "AliTPCExBBShape.h"
#include "TRandom.h"
#include "AliGeomManager.h"
#include "AliESDVertex.h"
#include "AliTPCcalibDB.h"
#endif


void AliTPCCorrectionDemo() {

  //
  // This is a Demo function of the general class AliTPCCorrection, which is used for 
  // general space point correction due to different effects.
  // The effects used in this Demo are:
  //   1. ExB twist - general offset of the TPC axis in comparison to the B field axis
  //   2. GG error (Gating Grid volt. error) - not perfectly aligned GG voltage (in terms of voltage)
  //   3. ExBBShape - B field shape correction of the secound order
  //
  // See class descriptions for further details 
  //
  // Authors: Magnus Mager, Stefan Rossegger, Jim Thomas
  //
  //
  // omegaTau (wt) of the langevin equation
  // This is a function of the drift vel., the magnetic and electric field
  // e.g. vd=2.6 cm/usc; Ez=400 V/cm; Bz=0.5 T
  // wt =  -10.0*(Bz*10)*vd/Ez = -0.325 

  Double_t vdrift = 2.6; // [cm/us]   // to be updated: per second (ideally)
  Double_t bzField = -0.5; // [Tesla] // to be updated: per run
  Double_t ezField = 400; // [V/cm]   // to be updated: never (hopefully)
  Double_t wt = -10.0 * (bzField*10) * vdrift / ezField ; 

  // Correction Terms for effective omegaTau; obtained by a laser calibration run
  Double_t T1 = 0.9;
  Double_t T2 = 1.5;

  AliMagF mag("mag","mag");

  AliTPCExBTwist twist;
  twist.SetXTwist(0.001);
  
  AliTPCGGVoltError GGerror;
  GGerror.SetDeltaVGGA(50.);
  GGerror.SetDeltaVGGC(50.);
  GGerror.InitGGVoltErrorDistortion();

  AliTPCExBBShape exb;
  exb.SetBField(&mag);

  TObjArray cs;
  cs.Add(&twist);
  cs.Add(&GGerror);
  cs.Add(&exb);

  AliTPCComposedCorrection cc;
  cc.SetCorrections(&cs);
  cc.SetOmegaTauT1T2(wt,T1,T2);
  //cc.SetMode(1);

  cc.Print("DA"); // Print used correction classes

  TCanvas *c=new TCanvas;  // Plots
  c->Divide(2,2);
  c->cd(1);twist.CreateHistoDRPhiinZR(1.,100,100)->Draw("surf2");
  c->cd(2);GGerror.CreateHistoDRPhiinZR(1.,100,100)->Draw("surf2");
  c->cd(3);exb.CreateHistoDRPhiinZR(1.)->Draw("surf2");
  c->cd(4);cc.CreateHistoDRPhiinZR(1.)->Draw("surf2");
}


void MakeDistortionMap(){
  //
  // make distortiona map example for specific transformation 
  //
  Int_t run=0;
  TTreeSRedirector *pcstream =  new TTreeSRedirector("distort.root");
  Double_t vdrift = 2.6; // [cm/us]   // to be updated: per second (ideally)
  Double_t bzField = -0.5; // [Tesla] // to be updated: per run
  Double_t ezField = 400; // [V/cm]   // to be updated: never (hopefully)
  Double_t wt = -10.0 * (bzField*10) * vdrift / ezField ; 

  // Correction Terms for effective omegaTau; obtained by a laser calibration run
  Double_t T1 = 0.9;
  Double_t T2 = 1.5;
  AliMagF* magF= new AliMagF("Maps","Maps", bzField/0.5, 1.);
  TGeoGlobalMagField::Instance()->SetField(magF);

  AliCDBManager::Instance()->SetDefaultStorage("local:////$ALICE_ROOT/OCDB");
  AliCDBManager::Instance()->SetRun(run);
  AliGeomManager::LoadGeometry();
  AliGeomManager::ApplyAlignObjsFromCDB("GRP ITS TPC");


  AliTPCExBTwist twist;
  twist.SetName("twistX001");
  twist.SetXTwist(0.001);
  
  AliTPCGGVoltError GGerror;
  GGerror.SetDeltaVGGA(50.);
  GGerror.SetDeltaVGGC(50.);
  GGerror.InitGGVoltErrorDistortion();
  GGerror.SetName("ggoffsetA50C50");
  AliTPCExBBShape exb;
  exb.SetBField(magF);
  exb.SetName("ExB");
  TObjArray cs;
  cs.Add(&twist);
  cs.Add(&GGerror);
  cs.Add(&exb);
  
  AliTPCComposedCorrection *cc = new AliTPCComposedCorrection;
  cc->SetCorrections(&cs);
  cc->SetOmegaTauT1T2(wt,T1,T2);
  //cc->SetMode(1);
  cc->SetName("composed");
  Double_t par[5]={0,0,0,0,0};
  Double_t cov[15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0};
  for (Double_t theta=-1.; theta<1; theta+=0.1){
    for (Double_t alpha=-3.14; alpha<3.14; alpha+=0.2){
      par[0]=0;
      par[1]=theta*85;
      par[2]=0;
      par[3]=theta;
      par[4]=(gRandom->Rndm()-0.5);
      AliExternalTrackParam trackIn(85,alpha,par,cov);
      //FitDistortedTrack(&trackIn, cc, 85, -1,pcstream);
      twist.FitDistortedTrack(trackIn, 85, -1,pcstream);
      GGerror.FitDistortedTrack(trackIn, 85, -1,pcstream);
      exb.FitDistortedTrack(trackIn, 85, -1,pcstream);
      cc->FitDistortedTrack(trackIn, 85, -1,pcstream);      
    }
  }
  delete pcstream;
}

void DrawDistortionMap(){
  //
  // Example -drawing of distortion maps
  //
  TFile f("distort.root");
  TTree *   fitDistorttwistX001= (TTree*)f.Get("fitDistorttwistX001");
  TTree *   fitDistortggoffsetA50C50=(TTree*)f.Get("fitDistortggoffsetA50C50");
  TTree *   fitDistortExB=(TTree*)f.Get("fitDistortExB");
  TTree *   fitDistortcomposed=(TTree*)f.Get("fitDistortcomposed");
  fitDistortExB->SetMarkerColor(1);
  fitDistorttwistX001->SetMarkerColor(2);
  fitDistortggoffsetA50C50->SetMarkerColor(4);
  fitDistortExB->SetMarkerStyle(20);
  fitDistorttwistX001->SetMarkerStyle(21);
  fitDistortggoffsetA50C50->SetMarkerStyle(22);
  //
  // example draw delta local y/r-phi as function of fi
  Int_t entries=0;

  entries = fitDistortExB->Draw("track1.fP[0]-track0.fP[0]:track0.fAlpha","","");
  TGraph * grY0= new TGraph(entries, fitDistortExB->GetV2(),fitDistortExB->GetV1());
  entries = fitDistortggoffsetA50C50->Draw("track1.fP[0]-track0.fP[0]:track0.fAlpha","","");
  TGraph * grY1= new TGraph(entries, fitDistortggoffsetA50C50->GetV2(),fitDistortggoffsetA50C50->GetV1());
  entries = fitDistorttwistX001->Draw("track1.fP[0]-track0.fP[0]:track0.fAlpha","","");
  TGraph * grY2= new TGraph(entries, fitDistorttwistX001->GetV2(),fitDistorttwistX001->GetV1());
  entries = fitDistortcomposed->Draw("track1.fP[0]-track0.fP[0]:track0.fAlpha","","");
  TGraph * grY3= new TGraph(entries, fitDistortcomposed->GetV2(),fitDistortcomposed->GetV1());
  grY0->SetMinimum(-0.4);
  grY0->SetMaximum(0.4);
  grY0->SetMarkerStyle(20), grY0->SetMarkerColor(1);
  grY1->SetMarkerStyle(21), grY1->SetMarkerColor(2);
  grY2->SetMarkerStyle(22), grY2->SetMarkerColor(4);
  grY3->SetMarkerStyle(24), grY3->SetMarkerColor(5); 
  grY0->Draw("ap");
  grY1->Draw("p");
  grY2->Draw("p");
  grY3->Draw("p");
}


void TestVertex(){
  //
  //
  //  .x ConfigCalibTrain.C(120829)

  AliTPCComposedCorrection * corrC = ( AliTPCComposedCorrection *)AliTPCcalibDB::Instance()->GetTPCComposedCorrection(0.5);
  AliTPCCorrection * corrT = (AliTPCCorrection *)corrC->GetCorrections()->FindObject("exb_twist");
  AliTPCCorrection * corrS = (AliTPCCorrection *)corrC->GetCorrections()->FindObject("ExB");
  
  TTreeSRedirector *pcstream = new TTreeSRedirector("vertexDistort.root");
  TTreeSRedirector *pcstreamS = new TTreeSRedirector("vertexDistortS.root");
  TTreeSRedirector *pcstreamT = new TTreeSRedirector("vertexDistortT.root");
  Double_t orgVertex[3] ;
  AliESDVertex aV,cV,aVO,cVO;
  for (Int_t iv=0; iv<100; iv++){
    printf("%d\n",iv);
    orgVertex[0]=gRandom->Gaus()*0.01;
    orgVertex[1]=gRandom->Gaus()*0.01;
    orgVertex[2]=gRandom->Gaus()*3;
    corrC->FastSimDistortedVertex(orgVertex,100, aV,aVO,cV,cVO,pcstream);
    corrS->FastSimDistortedVertex(orgVertex,100, aV,aVO,cV,cVO,pcstreamS);
    corrT->FastSimDistortedVertex(orgVertex,100, aV,aVO,cV,cVO,pcstreamT);
  }
  delete pcstream;
  delete pcstreamS;
  delete pcstreamT;
  TFile fC("vertexDistortC.root");
  TFile fT("vertexDistortT.root");
  TFile fS("vertexDistortS.root");
  TTree * treeT=(TTree*)fT.Get("vertex"); // twist distrotions
  TTree * treeS=(TTree*)fS.Get("vertex"); // shape distortions
  TCanvas *canvasD=new TCanvas("canvasD","canvasD");
  canvasD->Divide(2,2);
  canvasD->cd(1);
  treeT->SetLineColor(2);
  treeT->Draw("av.fPosition[0]-x>>hisATX(100,-0.2,0.2)","","");
  treeT->SetLineColor(4);
  treeT->Draw("cv.fPosition[0]-x>>hisCTX(100,-0.2,0.2)","","same");
  canvasD->cd(2);
  treeT->SetLineColor(2);
  treeT->Draw("av.fPosition[1]-y>>hisATY(100,-0.2,0.2)","","");
  treeT->SetLineColor(4);
  treeT->Draw("cv.fPosition[1]-y>>hisCTY(100,-0.2,0.2)","","same");
  //
  canvasD->cd(3);
  treeS->SetLineColor(2);
  treeS->Draw("av.fPosition[0]-x>>hisASX(100,-0.2,0.2)","","");
  treeS->SetLineColor(4);
  treeS->Draw("cv.fPosition[0]-x>>hisCSX(100,-0.2,0.2)","","same");
  canvasD->cd(4);
  treeS->SetLineColor(2);
  treeS->Draw("av.fPosition[1]-y>>hisASY(100,-0.2,0.2)","","");
  treeS->SetLineColor(4);
  treeS->Draw("cv.fPosition[1]-y>>hisCSY(100,-0.2,0.2)","","same");
  canvasD->SaveAs("vertexShift.ps");

}
Commit	Line	Data
f43f2c1e	1	/*
	2	gSystem->SetIncludePath("-I$ROOTSYS/include -I$ALICE_ROOT/ -I$ALICE_ROOT/include -I$ALICE_ROOT/STEER -I$ALICE_ROOT/TPC -I$ALICE_ROOT/ITS -I$ALICE_ROOT/TRD -I$ALICE_ROOT/TOF -I$ALICE_ROOT/RAW -I$ALICE_ROOT/STAT");
	3
	4	*/
	5	#if !defined(__CINT__) \|\| defined(__MAKECINT__)
	6	#include "THnSparse.h"
	7	#include "TLatex.h"
	8	#include "TCanvas.h"
	9	#include "TLegend.h"
	10	#include "TSystem.h"
	11	#include "TFile.h"
	12	#include "TChain.h"
	13	#include "TCut.h"
	14	#include "TH3.h"
	15	#include "TH2F.h"
	16	#include "TProfile3D.h"
	17	#include "TMath.h"
	18	#include "TVectorD.h"
	19	#include "TMatrixD.h"
	20	#include "TStatToolkit.h"
	21	#include "TTreeStream.h"
	22	#include "AliExternalTrackParam.h"
	23	#include "AliESDfriend.h"
	24	#include "AliTPCcalibTime.h"
	25	#include "TROOT.h"
	26	#include "AliXRDPROOFtoolkit.h"
	27	#include "AliTPCCorrection.h"
	28	#include "AliTPCExBTwist.h"
	29	#include "AliTPCGGVoltError.h"
	30	#include "AliTPCComposedCorrection.h"
	31	#include "AliTPCExBConical.h"
	32	#include "TPostScript.h"
	33	#include "TStyle.h"
	34	#include "AliTrackerBase.h"
	35	#include "TGraph.h"
	36	#include "AliCDBManager.h"
	37	#include "AliTPCExBBShape.h"
	38	#include "TRandom.h"
	39	#include "AliGeomManager.h"
ca58ed4e	40	#include "AliESDVertex.h"
ca58ed4e	41	#include "AliTPCcalibDB.h"
f43f2c1e	42	#endif
	43
	44
0116859c	45	void AliTPCCorrectionDemo() {
	46
	47	//
	48	// This is a Demo function of the general class AliTPCCorrection, which is used for
	49	// general space point correction due to different effects.
	50	// The effects used in this Demo are:
	51	// 1. ExB twist - general offset of the TPC axis in comparison to the B field axis
	52	// 2. GG error (Gating Grid volt. error) - not perfectly aligned GG voltage (in terms of voltage)
	53	// 3. ExBBShape - B field shape correction of the secound order
	54	//
	55	// See class descriptions for further details
	56	//
	57	// Authors: Magnus Mager, Stefan Rossegger, Jim Thomas
	58	//
	59	//
	60	// omegaTau (wt) of the langevin equation
	61	// This is a function of the drift vel., the magnetic and electric field
	62	// e.g. vd=2.6 cm/usc; Ez=400 V/cm; Bz=0.5 T
	63	// wt = -10.0(Bz10)*vd/Ez = -0.325
	64
	65	Double_t vdrift = 2.6; // [cm/us] // to be updated: per second (ideally)
	66	Double_t bzField = -0.5; // [Tesla] // to be updated: per run
	67	Double_t ezField = 400; // [V/cm] // to be updated: never (hopefully)
	68	Double_t wt = -10.0 * (bzField10) vdrift / ezField ;
	69
	70	// Correction Terms for effective omegaTau; obtained by a laser calibration run
	71	Double_t T1 = 0.9;
	72	Double_t T2 = 1.5;
	73
	74	AliMagF mag("mag","mag");
	75
	76	AliTPCExBTwist twist;
	77	twist.SetXTwist(0.001);
	78
	79	AliTPCGGVoltError GGerror;
	80	GGerror.SetDeltaVGGA(50.);
	81	GGerror.SetDeltaVGGC(50.);
	82	GGerror.InitGGVoltErrorDistortion();
	83
	84	AliTPCExBBShape exb;
	85	exb.SetBField(&mag);
	86
	87	TObjArray cs;
	88	cs.Add(&twist);
	89	cs.Add(&GGerror);
	90	cs.Add(&exb);
	91
	92	AliTPCComposedCorrection cc;
	93	cc.SetCorrections(&cs);
	94	cc.SetOmegaTauT1T2(wt,T1,T2);
f43f2c1e	95	//cc.SetMode(1);
0116859c	96
	97	cc.Print("DA"); // Print used correction classes
	98
	99	TCanvas *c=new TCanvas; // Plots
	100	c->Divide(2,2);
	101	c->cd(1);twist.CreateHistoDRPhiinZR(1.,100,100)->Draw("surf2");
	102	c->cd(2);GGerror.CreateHistoDRPhiinZR(1.,100,100)->Draw("surf2");
	103	c->cd(3);exb.CreateHistoDRPhiinZR(1.)->Draw("surf2");
	104	c->cd(4);cc.CreateHistoDRPhiinZR(1.)->Draw("surf2");
	105	}
e44441e2	106
	107
	108
	109	void MakeDistortionMap(){
	110	//
	111	// make distortiona map example for specific transformation
	112	//
	113	Int_t run=0;
	114	TTreeSRedirector *pcstream = new TTreeSRedirector("distort.root");
	115	Double_t vdrift = 2.6; // [cm/us] // to be updated: per second (ideally)
	116	Double_t bzField = -0.5; // [Tesla] // to be updated: per run
	117	Double_t ezField = 400; // [V/cm] // to be updated: never (hopefully)
	118	Double_t wt = -10.0 * (bzField10) vdrift / ezField ;
	119
	120	// Correction Terms for effective omegaTau; obtained by a laser calibration run
	121	Double_t T1 = 0.9;
	122	Double_t T2 = 1.5;
	123	AliMagF* magF= new AliMagF("Maps","Maps", bzField/0.5, 1.);
	124	TGeoGlobalMagField::Instance()->SetField(magF);
	125
	126	AliCDBManager::Instance()->SetDefaultStorage("local:////$ALICE_ROOT/OCDB");
	127	AliCDBManager::Instance()->SetRun(run);
	128	AliGeomManager::LoadGeometry();
	129	AliGeomManager::ApplyAlignObjsFromCDB("GRP ITS TPC");
	130
	131
	132	AliTPCExBTwist twist;
	133	twist.SetName("twistX001");
	134	twist.SetXTwist(0.001);
	135
	136	AliTPCGGVoltError GGerror;
	137	GGerror.SetDeltaVGGA(50.);
	138	GGerror.SetDeltaVGGC(50.);
	139	GGerror.InitGGVoltErrorDistortion();
	140	GGerror.SetName("ggoffsetA50C50");
	141	AliTPCExBBShape exb;
	142	exb.SetBField(magF);
	143	exb.SetName("ExB");
	144	TObjArray cs;
	145	cs.Add(&twist);
	146	cs.Add(&GGerror);
	147	cs.Add(&exb);
	148
	149	AliTPCComposedCorrection *cc = new AliTPCComposedCorrection;
	150	cc->SetCorrections(&cs);
	151	cc->SetOmegaTauT1T2(wt,T1,T2);
	152	//cc->SetMode(1);
	153	cc->SetName("composed");
	154	Double_t par[5]={0,0,0,0,0};
	155	Double_t cov[15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0};
	156	for (Double_t theta=-1.; theta<1; theta+=0.1){
	157	for (Double_t alpha=-3.14; alpha<3.14; alpha+=0.2){
	158	par[0]=0;
	159	par[1]=theta*85;
	160	par[2]=0;
	161	par[3]=theta;
	162	par[4]=(gRandom->Rndm()-0.5);
	163	AliExternalTrackParam trackIn(85,alpha,par,cov);
	164	//FitDistortedTrack(&trackIn, cc, 85, -1,pcstream);
f43f2c1e	165	twist.FitDistortedTrack(trackIn, 85, -1,pcstream);
	166	GGerror.FitDistortedTrack(trackIn, 85, -1,pcstream);
	167	exb.FitDistortedTrack(trackIn, 85, -1,pcstream);
	168	cc->FitDistortedTrack(trackIn, 85, -1,pcstream);
e44441e2	169	}
	170	}
	171	delete pcstream;
	172	}
	173
	174	void DrawDistortionMap(){
	175	//
	176	// Example -drawing of distortion maps
	177	//
	178	TFile f("distort.root");
	179	TTree * fitDistorttwistX001= (TTree*)f.Get("fitDistorttwistX001");
	180	TTree * fitDistortggoffsetA50C50=(TTree*)f.Get("fitDistortggoffsetA50C50");
	181	TTree * fitDistortExB=(TTree*)f.Get("fitDistortExB");
	182	TTree * fitDistortcomposed=(TTree*)f.Get("fitDistortcomposed");
	183	fitDistortExB->SetMarkerColor(1);
	184	fitDistorttwistX001->SetMarkerColor(2);
	185	fitDistortggoffsetA50C50->SetMarkerColor(4);
	186	fitDistortExB->SetMarkerStyle(20);
	187	fitDistorttwistX001->SetMarkerStyle(21);
	188	fitDistortggoffsetA50C50->SetMarkerStyle(22);
	189	//
	190	// example draw delta local y/r-phi as function of fi
	191	Int_t entries=0;
	192
	193	entries = fitDistortExB->Draw("track1.fP[0]-track0.fP[0]:track0.fAlpha","","");
	194	TGraph * grY0= new TGraph(entries, fitDistortExB->GetV2(),fitDistortExB->GetV1());
	195	entries = fitDistortggoffsetA50C50->Draw("track1.fP[0]-track0.fP[0]:track0.fAlpha","","");
	196	TGraph * grY1= new TGraph(entries, fitDistortggoffsetA50C50->GetV2(),fitDistortggoffsetA50C50->GetV1());
	197	entries = fitDistorttwistX001->Draw("track1.fP[0]-track0.fP[0]:track0.fAlpha","","");
	198	TGraph * grY2= new TGraph(entries, fitDistorttwistX001->GetV2(),fitDistorttwistX001->GetV1());
	199	entries = fitDistortcomposed->Draw("track1.fP[0]-track0.fP[0]:track0.fAlpha","","");
	200	TGraph * grY3= new TGraph(entries, fitDistortcomposed->GetV2(),fitDistortcomposed->GetV1());
	201	grY0->SetMinimum(-0.4);
	202	grY0->SetMaximum(0.4);
	203	grY0->SetMarkerStyle(20), grY0->SetMarkerColor(1);
	204	grY1->SetMarkerStyle(21), grY1->SetMarkerColor(2);
	205	grY2->SetMarkerStyle(22), grY2->SetMarkerColor(4);
	206	grY3->SetMarkerStyle(24), grY3->SetMarkerColor(5);
	207	grY0->Draw("ap");
	208	grY1->Draw("p");
	209	grY2->Draw("p");
	210	grY3->Draw("p");
	211	}
	212
ca58ed4e	213
	214	void TestVertex(){
	215	//
	216	//
	217	// .x ConfigCalibTrain.C(120829)
	218
	219	AliTPCComposedCorrection * corrC = ( AliTPCComposedCorrection *)AliTPCcalibDB::Instance()->GetTPCComposedCorrection(0.5);
	220	AliTPCCorrection * corrT = (AliTPCCorrection *)corrC->GetCorrections()->FindObject("exb_twist");
	221	AliTPCCorrection * corrS = (AliTPCCorrection *)corrC->GetCorrections()->FindObject("ExB");
	222
	223	TTreeSRedirector *pcstream = new TTreeSRedirector("vertexDistort.root");
	224	TTreeSRedirector *pcstreamS = new TTreeSRedirector("vertexDistortS.root");
	225	TTreeSRedirector *pcstreamT = new TTreeSRedirector("vertexDistortT.root");
	226	Double_t orgVertex[3] ;
	227	AliESDVertex aV,cV,aVO,cVO;
	228	for (Int_t iv=0; iv<100; iv++){
	229	printf("%d\n",iv);
	230	orgVertex[0]=gRandom->Gaus()*0.01;
	231	orgVertex[1]=gRandom->Gaus()*0.01;
	232	orgVertex[2]=gRandom->Gaus()*3;
	233	corrC->FastSimDistortedVertex(orgVertex,100, aV,aVO,cV,cVO,pcstream);
	234	corrS->FastSimDistortedVertex(orgVertex,100, aV,aVO,cV,cVO,pcstreamS);
	235	corrT->FastSimDistortedVertex(orgVertex,100, aV,aVO,cV,cVO,pcstreamT);
	236	}
	237	delete pcstream;
	238	delete pcstreamS;
	239	delete pcstreamT;
	240	TFile fC("vertexDistortC.root");
	241	TFile fT("vertexDistortT.root");
	242	TFile fS("vertexDistortS.root");
	243	TTree * treeT=(TTree*)fT.Get("vertex"); // twist distrotions
	244	TTree * treeS=(TTree*)fS.Get("vertex"); // shape distortions
	245	TCanvas *canvasD=new TCanvas("canvasD","canvasD");
	246	canvasD->Divide(2,2);
	247	canvasD->cd(1);
	248	treeT->SetLineColor(2);
	249	treeT->Draw("av.fPosition[0]-x>>hisATX(100,-0.2,0.2)","","");
	250	treeT->SetLineColor(4);
	251	treeT->Draw("cv.fPosition[0]-x>>hisCTX(100,-0.2,0.2)","","same");
	252	canvasD->cd(2);
	253	treeT->SetLineColor(2);
	254	treeT->Draw("av.fPosition[1]-y>>hisATY(100,-0.2,0.2)","","");
	255	treeT->SetLineColor(4);
	256	treeT->Draw("cv.fPosition[1]-y>>hisCTY(100,-0.2,0.2)","","same");
	257	//
	258	canvasD->cd(3);
	259	treeS->SetLineColor(2);
	260	treeS->Draw("av.fPosition[0]-x>>hisASX(100,-0.2,0.2)","","");
	261	treeS->SetLineColor(4);
	262	treeS->Draw("cv.fPosition[0]-x>>hisCSX(100,-0.2,0.2)","","same");
	263	canvasD->cd(4);
	264	treeS->SetLineColor(2);
	265	treeS->Draw("av.fPosition[1]-y>>hisASY(100,-0.2,0.2)","","");
	266	treeS->SetLineColor(4);
	267	treeS->Draw("cv.fPosition[1]-y>>hisCSY(100,-0.2,0.2)","","same");
	268	canvasD->SaveAs("vertexShift.ps");
	269
	270	}