// 24 luglio 2007
//----------------------------
//	anna.parravicini@cnao.it
//----------------------------

/*
routine per plottare distribuzione particelle
nello spazio fasi ed ellisse di Twiss

>>>> Difference with previous versions:

*the reference ellipse is plotted together with the RMS ellipse
*geometrical emittance plot added
*TProfile instead of TH2F, in order to have the average z-entry
 rather than the sum of z-entries.
*x-bins [mm] are centred on the slit median axis width. namely,
 the bin is drawn from xo-(slit step)/2 to xo+(slit step)/2
*main parameters written on the graphs
*wrt v4, aggiungo griglia ai pad e banner per output
*wrt v5, bunner spostato a sinistra e scritte molto piu' cose,
  plottati dati prima del centraggio
*wrt v6-7, aggiungo isocurva relative a Ithreshold sul grafico emitt geom.
 scrivo nel bunner la percentuale di corrente di fascio contenuta nelle
 ellisse rms e reference
*wrt v8, modificato per lavorare con emittGeom_v1	

>>>> To be done:

* aggiungere logo CNAO
* si potrebbero cercare i valori estremi presenti negli array con corrente nn nulla
 e settare i range dei grafici in modo automatico
* decidere come settare il valore dell'errore sull'emittanza geometrica

//----------------------------------------------------------*/

#include "TGraph.h"
#include "TStyle.h"
#include "TMath.h"
#include "TCanvas.h"
#include "TTree.h"
#include "TProfile2D.h" //user's guide, page 44. See prova_TH2F_TProfile2D.C per esempio.

int Main(){

	// === VARIABLES DECLARATION ===================================================
	Int_t Nentries; //events number per branch

	Float_t pos[100000];//[mm] - array da riempire con Branches
	Float_t div[100000];//[mrad] - array da riempire con Branches
	Float_t charge[100000];//[arbitrary units] - array da riempire con Branches

	Float_t pos_t=0;//variabili d'appoggio per leggere Tree
    Float_t div_t=0;
	Float_t charge_t=0;

	Double_t Xmax;//extreme values for x and xp arrays
	Double_t Xmin;
	Double_t XPmax;
	Double_t XPmin;
	Double_t Imax;
	Double_t Imin;

	//RMS ellipse parameters
    Double_t alpha;
    Double_t beta;
    Double_t gamma;
    Double_t emi;
    Double_t limit;

	Int_t n=44;
	Double_t ell[88];
    Double_t x[88];
    Double_t step;

	//particles barycenter
	Double_t x_center;//[mm]
	Double_t xp_center;//[mrad]
	Int_t center; //0=data not centred; 1=data centred

	//bias and noise thresholds
	Int_t ThrB; // % of max entry current
	Int_t ThrZ; // % of max entry current

	//reference ellipse parameters
	Double_t alpha_ref;
    Double_t beta_ref;
    Double_t gamma_ref;
    Double_t emi_ref;
    Double_t limit_ref;

	Double_t ell_ref[88];
    Double_t x_ref[88];
    Double_t step_ref;

	//inside beam percentage
	Float_t rms_IN; //beam charge percentage inside RMS ellipse
	Float_t ref_IN;//beam charge percentage inside reference ellipse

	//geometrical ellipse threshold
	Double_t Ithreshold;
	Bool_t overloop; //in case it doesn't reach the required accuracy
	Float_t P; //we want compute 90% emittance
	Float_t Peff; //effective percentage of the emittance computation
	Double_t emiP;//[mm*mrad] geometrical emittance output value
	Double_t emiRes;//[mm*mrad] geometrical emittance error
	Double_t maxDiscr;//[mm*mrad] maximum discrepancy among required (P) and computed (Peff) emittance

	//scanning parameters
	Int_t sl_Nsteps;//number of slit steps
	Float_t sl_step;//slit step amplitude
	Int_t ws_points;// ws samples per profile, at constant velocity

	//measurement name
	string name;
	string outname;
	string bunner;

	//general settings
	gStyle->SetOptStat(kFALSE);//not to plot the statistic box
	gStyle->SetPalette(1);
	//gStyle->SetOptDate(31);//rimane nascosta sotto ai pad, non si vede

	// === CREO UN CANVAS in cui fare i plot e lo divido in 3 =======================

	TCanvas *c=new TCanvas ("c", "Phase Space Plots",1200, 800);
	c->Divide(2,2);

	// ==== CARICO PUNTI SP FASI DA FILE (Il file deve stare in C:\root) =============

	// --- from file to the tree-------------------------------

	// per fare histo 2D da file con 3 colonne, di cui la terza e' il "peso"
	// da dare alla cella
	TTree *t=new TTree ("t", "albero per caricare file di 3 colonne");
	t->ReadFile("meas_distribution.txt", "x:xp:I");

	Nentries=(Int_t)t->GetEntries(); //numero di eventi per ramo
	Xmax=t->GetMaximum("x");
	Xmin=t->GetMinimum("x");
	XPmax=t->GetMaximum ("xp");
	XPmin=t->GetMinimum ("xp");
	Imax=t->GetMaximum ("I");
	Imin=t->GetMinimum ("I");

	// --- from tree to the arrays------------------------------

	//se mi bastano histo e non voglio grafico,
	//non e' necessario leggere l'albero e riempire gli array
	
	//before reading branches, we must give them an address
	t->SetBranchAddress("x",&pos_t);
	t->SetBranchAddress("xp",&div_t);
	t->SetBranchAddress("I",&charge_t);
	
	//fill the array from the branches
	for (Int_t j=0;j<Nentries;j++){
		t->GetEntry(j);
		pos[j]=pos_t;
		div[j]=div_t;
		charge[j]=charge_t;
	}

	// === CARICO PARAMETRI DI TWISS e preparo i punti per disegnare l'ellisse ==========

	ifstream infile("twisspara.txt");
    Double_t param[24];
	// open the twisspara.txt file to understand the list of parameters written there

    for (Int_t i=0;i<22;i++)
	   infile>>param[i];
	infile>>name;
	name=name.substr(13);
	infile>>center;
	
	//RMS ellipse parameters
    alpha=param[0];
    beta=param[1];
    gamma=(1+alpha*alpha)/beta;
    emi=param[2];
    limit=sqrt(emi*beta);

	//reference ellipse parameters
	alpha_ref=param[3];
    beta_ref=param[4];
    gamma_ref=(1+alpha_ref*alpha_ref)/beta_ref;
    emi_ref=param[5];
    limit_ref=sqrt(emi_ref*beta_ref);

	//particles barycenter
	x_center=param[6];
	xp_center=param[7];

	//bias and noise thresolds
	Double_t thr;
	thr=param[18];
	ThrB= (int) thr*100; 
	thr=param[19];
	ThrZ= (int) thr*100;

	//beam percentage inside ellipses
	rms_IN=param[20];//beam current percentage inside rms ellipse
	ref_IN=param[21];//insidea reference ellipse

	//scanning settings
	sl_Nsteps=param[9];
	sl_step=param[10];
	ws_points=param[11];

	//geometrical ellipse
	Ithreshold=param[8];
	
	overloop=param[12];
	if (overloop)
		cout<<"Geom emitt is far from user's required accuracy."<<endl; //da aggiungere anche sul plot
	
	P=param[13]*100;// %
	Peff=param[14]*100;// %
	emiP=param[15];//[mm*mrad]
	emiRes=param[16];//[mm*mrad]
	maxDiscr=param[17];//[mm*mrad]

	//writing the RMS ellipse 
	step=limit/(n/2);

    for (Int_t i=1;i<n;i++){ //riempie gli array con i punti dell' ellisse
		//in the border points, it could happen the number under sqrt is not zero but slightly
		//negative. I put the number unders sqrt 
	   x[i]=-limit+step*i;
	   ell[i]=(-2*alpha*x[i]+sqrt(fabs((2*alpha*x[i])*(2*alpha*x[i])-4*beta*(gamma*x[i]*x[i]-emi))))/(2*beta);
 	   x[i+n]=limit-step*i;
	   ell[i+n]=(-2*alpha*x[i+n]-sqrt(fabs((2*alpha*x[i+n])*(2*alpha*x[i+n])-4*beta*(gamma*x[i+n]*x[i+n]-emi))))/(2*beta);
    }
	//riempio separatamente i punti estremi perche' per questi puo' capitare che il radicando sia molto
	//piccolo, ma negativo, dando quindi soluzioni non reali
   	x[0]=-limit;
	ell[0]=-2*alpha*x[0]/(2*beta);
	x[n]=limit;
	ell[n]=-2*alpha*x[n]/(2*beta);

	//writing the reference ellipse
	step_ref=limit_ref/(n/2);
	//cout<<"limit_ref="<<limit_ref<<endl;
   
    for (Int_t i=0;i<n;i++){ //riempie gli array con i punti dell' ellisse di riferimento
	   x_ref[i]=-limit_ref+step_ref*i;
	  //cout<<"x_ref[1]="<<x_ref[1]<<endl;
	   ell_ref[i]=(-2*alpha_ref*x_ref[i]+sqrt(fabs((2*alpha_ref*x_ref[i])*(2*alpha_ref*x_ref[i])-4*beta_ref*(gamma_ref*x_ref[i]*x_ref[i]-emi_ref))))/(2*beta_ref);
 	  // cout<<"ell_ref["<<i<<"]="<<ell_ref[i]<<endl;
	   x_ref[i+n]=limit_ref-step_ref*i;
	   ell_ref[i+n]=(-2*alpha_ref*x_ref[i+n]-sqrt(fabs((2*alpha_ref*x_ref[i+n])*(2*alpha_ref*x_ref[i+n])-4*beta_ref*(gamma_ref*x_ref[i+n]*x_ref[i+n]-emi_ref))))/(2*beta_ref);
    }
   	
	// === costruiamo i GRAFICI che RIPRODUCANO l'ELLISSE ===============================
	
	TGraph *gEll=new TGraph(2*n,x,ell);
	gEll->SetLineColor(1);
    //gEll->SetLineStyle(2);//2=dash, 4=dash-dot
	gEll->SetLineWidth(2);
	
	TGraph *gEllRef=new TGraph(2*n,x_ref,ell_ref);
	gEllRef->SetLineColor(1);
	gEllRef->SetLineStyle(2);//2=dash, 4=dash-dot
	gEllRef->SetLineWidth(2);

	// === PHASE SPACE DISTRIBUTION (x, xp, I) PLOT ===============================================================
	c->cd(1);
	gPad->SetPad(0.02,0.5,0.32,1);
	
	TProfile2D *h3D = new TProfile2D("h3d","#font[42]{Phase space distribution}", sl_Nsteps,Xmin-sl_step/2, Xmax+sl_step/2,ws_points, XPmin, XPmax+0.0001*XPmax);
	//TH2F *h3d = new TH2F("h3d","#font[42]{Phase space distribution}", sl_steps,Xmin, Xmax,ws_points, XPmin, XPmax);
	
	//TH2F plots the sum of all the z-entries of each bin;
	//TProfile2D plots the average of all the z-entries of each bin.

	for (int j=0;j<Nentries;j++)
		h3d->Fill(pos[j],div[j],charge[j]);
//	h3d->GetYaxis()->SetRange(ws_points/3,ws_points*2/3); //range di bin da plottare
	h3d->GetXaxis()->SetTitle("#font[42]{x[mm]}");
	h3d->GetYaxis()->SetTitle("#font[42]{xp[mrad]}");
	h3d->Draw("lego2z fb bb");
	
	//fb=remove front surface; bb=remove back surface of the cartesian axis

	// === NOT CENTRED DATA ==================================================
	c->cd(2);
	gPad->SetPad(0.34,0.5,0.64,1);
	gPad->SetGrid(); //gpad e' il puntatore al pad attivo

	TProfile2D *hnocenter = new TProfile2D("hnocenter","#font[42]{Phase space distribution (before centring data)}", sl_Nsteps,Xmin-sl_step/2+x_center, Xmax+sl_step/2+x_center,ws_points, XPmin+xp_center, XPmax+0.001*XPmax+xp_center);
	for (int j=0;j<Nentries;j++)
		hnocenter->Fill(pos[j]+x_center,div[j]+xp_center,charge[j]);
//	hnocenter->GetYaxis()->SetRange(ws_points/3,ws_points*2/3); //range di bin da plottare
	hnocenter->GetXaxis()->SetTitle("#font[42]{x[mm]}");
	hnocenter->GetYaxis()->SetTitle("#font[42]{xp[mrad]}");
	gStyle->SetOptStat(kFALSE);//not to plot the statistic box
	//hnocenter->Draw("cont1z"); //questo opzione riporta le isolinee, senza etichetta di valore	
	hnocenter->Draw("colz");

	/*** SCATTER PLOT (not very useful at now)
	// SCATTER PLOT 3D (non usa i colori, quindi non si capisce molto)
	TH3F *h3 = new TH3F("h3","#font[42]{Phase space scatter plot}", sl_Nsteps,Xmin, Xmax,ws_points, XPmin, XPmax, 100, Imin, Imax);
	for (int j=0;j<Nentries;j++)
		h3->Fill(pos[j],div[j],charge[j]);
	h3->GetYaxis()->SetRange(ws_points/3,ws_points*2/3); //range di bin da plottare
	h3->GetXaxis()->SetTitle("#font[42]{x[mm]}");
	h3->GetYaxis()->SetTitle("#font[42]{xp[mrad]}");
	h3->GetZaxis()->SetTitle("#font[42]{I[a.u.]}");
	h3->SetMarkerStyle(20);//pag 111 user's guide
    h3->SetMarkerSize(0.5);//pag 111 user's guide
	h3->Draw("");
	*/

	// === PHASE SPACE TOP-VIEW WITH RMS AND REFERENCE ELLIPSES  ==================================================
	c->cd(3);
	gPad->SetPad(0.02,0,0.32,0.5);
	gPad->SetGrid(); //gpad e' il puntatore al pad attivo
	//gPad->SetFrameFillColor(33);

	// histo 2D, with weigth I for each point.
	// If more points have the same x and xp, their I-values are sum up; it 
	// follows that a I value larger than Imax is displayed: it's due to the binning!
	TProfile2D *h = new TProfile2D("h","#font[42]{Statistical emittance and ref ellipse}", sl_Nsteps,Xmin-sl_step/2, Xmax+sl_step/2,ws_points, XPmin, XPmax+0.001*XPmax);
	for (int j=0;j<Nentries;j++)
		h->Fill(pos[j],div[j],charge[j]);
//	h->GetYaxis()->SetRange(ws_points/3,ws_points*2/3); //range di bin da plottare
	h->GetXaxis()->SetTitle("#font[42]{x[mm]}");
	h->GetYaxis()->SetTitle("#font[42]{xp[mrad]}");
	gStyle->SetOptStat(kFALSE);//not to plot the statistic box
	h->Draw("colz");

	// add the legend
	leg = new TLegend(0.47,0.11,0.9,0.25);
	leg->AddEntry(gEll,"RMS ellipse","l");//"l" use the line that markes the gEll graph
	leg->AddEntry(gEllRef,"Reference ellipse","l");
	leg->Draw();

	//Add statistical emitt value on the plot
	TPaveText *ptrms=new TPaveText (0.15,0.81,0.9,0.91,"NDC");
	//NDC: normalized coordinates (x1,y1,x2,y2)
	char text[100];
	sprintf (text,"#font[42]{#splitline{alpha RMS= %.3g - beta RMS= %.3g m}{Emitt RMS= %.3g #pi#upointmm#upointmrad}}", alpha, beta,emi);
	TText *trms=ptrms->AddText(text);
	trms->SetTextColor(1);
	trms->SetTextSize(0.04);
	trms->SetTextAlign(12);//=H and V centred

	gEll->Draw("C");//non riportiamo gli assi ("A") perche' gia' presenti. 
	gEllRef->Draw("C");
	ptrms->Draw();

	// === GEOMETRICAL EMITTANCE PLOT ============================================================================
	c->cd(4);
	gPad->SetPad(0.34,0,0.64,0.5);
	gPad->SetGrid();

	// param[8]=Ithreshold, namely, all the bins with current >Ithreshold
	// contribute to the geometrical emittance, at the required percentage.
	// I would like a top-view of the 3D phase-space plot, with a different fill color
	// for bins over and under threshold	
	// Only pints with I >Ithreshold are taken into account

	char hgeom_name[55];
	sprintf(hgeom_name,"#font[42]{Geometrical emittance surface @ %.0f%%}",P);

	TProfile2D *hgeom = new TProfile2D("hgeom",hgeom_name,sl_Nsteps,Xmin-sl_step/2, Xmax+sl_step/2,ws_points, XPmin, XPmax+0.001*XPmax);
	for (int j=0;j<Nentries;j++){
		if (charge[j]>Ithreshold)
			hgeom->Fill(pos[j],div[j],charge[j]);
	}

	//Add geom emitt value on the plot
	TPaveText *ptgeom=new TPaveText (0.25,0.81,0.9,0.91,"NDC"); //NDC: normalized coordinates (x1,y1,x2,y2)
	char text[100];
	sprintf (text,"#font[42]{#splitline{Geom Emitt=(%.2f#pm%.2f)#pi#upointmm#upointmrad}{(effective percentage: %.2f%%)}}",emiP,emiRes,Peff);
	TText *tgeom=ptgeom->AddText(text);
	tgeom->SetTextColor(1);
	tgeom->SetTextSize(0.04);
	tgeom->SetTextAlign(12);//=H and V centred

	//Add warning in case of bad geom emitt computation
	TPaveText *ptwarn=new TPaveText (0.15,0.15,0.7,0.3,"NDC");	//NDC: normalized coordinates (x1,y1,x2,y2)
	ptwarn->SetBorderSize(1);
	char text[100];
	sprintf (text,"#font[102]{#splitline{Geom emitt accuracy far}{from user request}}");
	TText *twarn=ptwarn->AddText(text);
	twarn->SetTextColor(2);
	twarn->SetTextSize(0.04);
	twarn->SetTextAlign(12);//=H and V centred

	//Add contour relative to Ithreshold
	//define contours value
	Double_t contours[1];
	contours[0]=Ithreshold;

	//firstly, we must plot the histo to extract the contour
	hgeom->SetContour(1, contours);//se funziona, mettere qui direttamente Ithreshold
	hgeom->Draw("cont z list"); //questo disegna le isolinee, non le numera pero'
	c->Update();//e' necessario plottare il grafico per estrarre le line di livello

	//get contours
	TObjArray *cont_number = (TObjArray*)gROOT->GetListOfSpecials()->FindObject("contours");
	TList *contLevel = NULL;
	TGraph * curv = NULL;
	Int_t nGraphs =0;
	Int_t TotalConts =0;

	if (cont_number==NULL){
		printf ("no contours extracted \n");
		TotalConts=0;
		return;
	}
	else {
		TotalConts=cont_number->GetSize();
		//TotalConts e' il numero di "contorni/isolinee" che trova nell'ultimo histo disegnato
	}
	printf ("Total contours =%d \n", TotalConts);

	for (int j=0; j<TotalConts; j++){

		contLevel=(TList*)cont_number->At(j); 
		//contLevel e' la lista di grafici (linee disgiunte) che formano il j-th contorno (isolinea)

		printf("Contour %d has %d Graphs\n", j, contLevel->GetSize());
		// contLevel->GetSize() dice il numero di grafici (linee disgiunte) che formano il j-th contorno
	}
	nGraphs=0;

	//TCanvas *cprova =new TCanvas ("c1", "Contour List", 610,0,600,600);
	TProfile2D *hgeom2 = new TProfile2D("hgeom",hgeom_name,sl_Nsteps,Xmin-sl_step/2, Xmax+sl_step/2,ws_points, XPmin, XPmax+0.001*XPmax);
	for (int j=0;j<Nentries;j++){
		if (charge[j]>Ithreshold)
			hgeom2->Fill(pos[j],div[j],charge[j]);
	}		
	hgeom2->Draw("colz"); //disegna il grafico su cui saranno plottate le isolinee
	//	hgeom2->GetYaxis()->SetRange(ws_points/3,ws_points*2/3); //range di bin da plottare
	hgeom2->GetXaxis()->SetTitle("#font[42]{x[mm]}");
	hgeom2->GetYaxis()->SetTitle("#font[42]{xp[mrad]}");

	//write isocurves labels
	TLatex la;
	la.SetTextSize(0.035);
	la.SetTextFont(42);
	char val[20];
	Double_t xo,yo,zo;

	for (int j=0; j<TotalConts; j++){ //TotalCont e' il numero totale di isolinee/contorni
		contLevel= (TList*)cont_number->At(j); //contLevel e' la lista dei grafici/curve dell'isolinea jth

		zo=contours[j];
		printf ("Isocurve passing by the contour Z=%f\n", zo);
		cout<<"has "<<contLevel->GetSize()<<" graphs"<<endl;
		
		curv=(TGraph*)contLevel->First();
		//scrive in "curv" il primo grafico (linea disgiunta) che forma il contorno

		for (int i=0;i<contLevel->GetSize();i++){//cicla fino al numero di grafici che formano l'isolinea

			curv->GetPoint(0,xo,yo);

			nGraphs++;
			
			printf ("\t Graph: %d -- %d elements \n", nGraphs,curv->GetN());
			//cosa fa GetN? Conta le entrate che stanno sulla linea "curv"?

			curv->SetLineWidth(2.5);
			curv->Draw("C");

			sprintf(val, "%.3g", zo);
			
			//la.DrawLatex(xo,yo,val);
			// Here above line plots isocurve value on the plot itself. Most of the time
			// it creates confusion with no help...so I keep the line commented at now
			
			curv=(TGraph*)contLevel->After(curv);// prende il grafico successivo della TList
		}//end i-loop

	}//end j-loop
	
	// add the legend
	leg2 = new TLegend(0.40,0.11,0.9,0.20);
	leg2->SetTextSize(0.035);
	leg2->AddEntry(curv,"Current threshold contour","l");//"l" use the line that markes the gEll graph
	leg2->Draw();

	c->Update();

	ptgeom->Draw();
	if (overloop)
		ptwarn->Draw();

	// === OUTPUT BANNERS ==================================================
	// Add to the canvas the TPave where writing the measurement information
	c->cd();

	//let's number the pads
	TPaveText *pad1=new TPaveText (0,0.5,0.019,1,"NDC"); //NDC: normalized coordinates (x1,y1,x2,y2)
	pad1->SetBorderSize(1);
	pad1->SetFillColor(29);
	TText *tpad1=pad1->AddText("PLOT 1");
	tpad1->SetTextColor(1);
	tpad1->SetTextSize(0.025);//pixel size= textSize*canvasWidth 
	tpad1->SetTextFont(102);
	tpad1->SetTextAlign(22);
	tpad1->SetTextAngle(90);
	pad1->Draw();

	TPaveText *pad2=new TPaveText (0.32,0.5,0.339,1,"NDC");
	pad2->SetBorderSize(1);
	pad2->SetFillColor(29);
	TText *tpad2=pad2->AddText("PLOT 2");
	tpad2->SetTextColor(1);
	tpad2->SetTextSize(0.025);//pixel size= textSize*canvasWidth 
	tpad2->SetTextFont(102);
	tpad2->SetTextAlign(22);
	tpad2->SetTextAngle(90);
	pad2->Draw();

	TPaveText *pad3=new TPaveText (0,0,0.019,0.5,"NDC");
	pad3->SetBorderSize(1);
	pad3->SetFillColor(29);
	TText *tpad3=pad3->AddText("PLOT 3");
	tpad3->SetTextColor(1);
	tpad3->SetTextSize(0.025);//pixel size= textSize*canvasWidth 
	tpad3->SetTextFont(102);
	tpad3->SetTextAlign(22);
	tpad3->SetTextAngle(90);
	pad3->Draw();

	TPaveText *pad4=new TPaveText (0.32,0,0.339,0.5,"NDC");
	pad4->SetBorderSize(1);
	pad4->SetFillColor(29);
	TText *tpad4=pad4->AddText("PLOT 4");
	tpad4->SetTextColor(1);
	tpad4->SetTextSize(0.025);//pixel size= textSize*canvasWidth 
	tpad4->SetTextFont(102);
	tpad4->SetTextAlign(22);
	tpad4->SetTextAngle(90);
	pad4->Draw();

	//..................

	//let's write the output information
	TPaveText *ptbunner=new TPaveText (0.64,0,1,1,"NDC");
	ptbunner->SetBorderSize(1);
	ptbunner->SetFillColor(29);

	TText *t1=ptbunner->AddText("EMITTANCE MEASUREMENT PROGRAM");
	t1->SetTextColor(1);
	t1->SetTextSize(0.03);//pixel size= textSize*canvasWidth 
	t1->SetTextFont(102);
	t1->SetTextAlign(22);

	TText *t2=ptbunner->AddText("");
	t2->SetTextColor(1);
	t2->SetTextSize(0.015);
	t2->SetTextFont(102);
	t2->SetTextAlign(23);

	bunner="#splitline{Input file name:}{";
	bunner.append(name);
	bunner.append ("}");
	TText *t3=ptbunner->AddText(bunner.c_str());
	t3->SetTextColor(1);
	t3->SetTextSize(0.02); 
	t3->SetTextFont(102);
	t3->SetTextAlign(13);

	char elab[200];
	if (center)
		sprintf (elab,"#splitline{Data elaboration parameters:}{Bias thr %d%%, Noise thr %d%%, centred}", ThrB, ThrZ);
	else 
		sprintf (elab,"#splitline{Data elaboration parameters:}{Bias thr %d%%, Noise thr %d%%, not-centred}", ThrB, ThrZ);
	TText *t4=ptbunner->AddText(elab);
	t4->SetTextColor(1);
	t4->SetTextSize(0.02);//pixel size= textSize*canvasWidth 
	t4->SetTextFont(102);
	t4->SetTextAlign(13);

	TText *t10=ptbunner->AddText("#diamond PLOT 1");
	t10->SetTextColor(1);
	t10->SetTextSize(0.025);
	t10->SetTextFont(102);
	t10->SetTextAlign(13);
	char description101[200];
	sprintf (description101,"#splitline{3D plot of elaborated data used to compute}{RMS and Geom emittance.}");
	TText *t101=ptbunner->AddText(description101);
	t101->SetTextColor(1);
	t101->SetTextSize(0.02); 
	t101->SetTextFont(102);
	t101->SetTextAlign(13);

	TText *t20=ptbunner->AddText("#diamond PLOT 2");
	t20->SetTextColor(1);
	t20->SetTextSize(0.025);
	t20->SetTextFont(102);
	t20->SetTextAlign(13);			
	char description201[200];
	sprintf (description201,"#splitline{Top-view of not centred data.}{Barycenter: x=%.3f mm, xp=%.3f mrad}",x_center, xp_center);
	TText *t201=ptbunner->AddText(description201);
	t201->SetTextColor(1);
	t201->SetTextSize(0.02);
	t201->SetTextFont(102);
	t201->SetTextAlign(13);

	TText *t30=ptbunner->AddText("#diamond PLOT 3");
	t30->SetTextColor(1);
	t30->SetTextSize(0.025);
	t30->SetTextFont(102);
	t30->SetTextAlign(13);			
	char description301[200];
	sprintf (description301,"#splitline{Top-view of elaborated data. RMS and}{reference ellipses plotted on it.}");
	TText *t301=ptbunner->AddText(description301);
	t301->SetTextColor(1);
	t301->SetTextSize(0.02);
	t301->SetTextFont(102);
	t301->SetTextAlign(13);
	char description302[200];
	sprintf (description302,"#splitline{Reference ellipse: beam charge inside %d%%,}{alpha=%.3g, beta=%.3g m, emitt=%.3g #pi#upointmm#upointmrad}", ref_IN,alpha_ref, beta_ref, emi_ref);
	TText *t302=ptbunner->AddText(description302);
	t302->SetTextColor(1);
	t302->SetTextSize(0.02);
	t302->SetTextFont(102);
	t302->SetTextAlign(13);
	char description303[200];
	sprintf (description303,"#splitline{RMS ellipse: beam charge inside %d%%, }{alpha=%.3g, beta=%.3g m, emitt=%.3g #pi#upointmm#upointmrad}",rms_IN, alpha, beta, emi);
	TText *t303=ptbunner->AddText(description303);
	t303->SetTextColor(1);
	t303->SetTextSize(0.02);
	t303->SetTextFont(102);
	t303->SetTextAlign(13);

	TText *t40=ptbunner->AddText("#diamond PLOT 4");
	t40->SetTextColor(1);
	t40->SetTextSize(0.025);
	t40->SetTextFont(102);
	t40->SetTextAlign(13);			
	char description401[200];
	sprintf (description401,"#splitline{Top-view of elaborated data contributing}{to the Geom emitt computation.}");
	TText *t401=ptbunner->AddText(description401);
	t401->SetTextColor(1);
	t401->SetTextSize(0.02); 
	t401->SetTextFont(102);
	t401->SetTextAlign(13);
	char description402[200];
	sprintf (description402,"#splitline{Required=%.1f%%, Effective=%.1f%%}{Current threshold=%.3g uA}", P,Peff,Ithreshold);
	TText *t402=ptbunner->AddText(description402);
	t402->SetTextColor(1);
	t402->SetTextSize(0.02); 
	t402->SetTextFont(102);
	t402->SetTextAlign(13);
	char description403[200];
	sprintf (description403,"#splitline{Maximum discrepancy between %.1f%% and %.1f%%}{geom emitt is %.2f#pi#upointmm#upointmrad}", P,Peff,maxDiscr);
	TText *t403=ptbunner->AddText(description403);
	t403->SetTextColor(1);
	t403->SetTextSize(0.02); 
	t403->SetTextFont(102);
	t403->SetTextAlign(13);
	char description404[200];
	sprintf (description404,"Geom Emitt at %.1f%%=(%.2f#pm%.2f)#pi#upointmm#upointmrad", Peff,emiP,emiRes);
	TText *t404=ptbunner->AddText(description404);
	t404->SetTextColor(1);
	t404->SetTextSize(0.02); 
	t404->SetTextFont(102);
	t404->SetTextAlign(13);

	ptbunner->Draw();

	// === SAVE THE OUTPUT ===========================================================================
	c->cd();
	
	outname="C:\\measOutput\\";
	outname.append(name);   
	outname.append("_plot.gif");
	c->SaveAs(outname.c_str());
	
	return 0;
}