// compilation:
// g++ `root-config --cflags` sandbox.cpp `root-config --glibs` -std=c++11
// run:
// /a.out

#include <TCanvas.h>
#include <TF1.h>
#include <TH1.h>
#include <TH2.h>
#include <TRandom.h>
#include <TRandom3.h>
#include <TFile.h>
#include <cmath>
#include <TVirtualFFT.h>
#include <TApplication.h>
#include <TLorentzVector.h>
#include <iostream>
#include <TGeoMatrix.h>
#include <map>
#include <cctype>
#include <utility>
#include <sstream>
#include <TGraphErrors.h>
#include <TGraphAsymmErrors.h>
#include <TTree.h>
#include <TMarker.h>
#include <TLine.h>
#include <TAxis.h>
#include <vector>


using namespace std;

const double t_min = 1e-6;
const double t_max = 1e-3;

const double energy = 1.5;

const double mass = 0.938;

const int n_steps = 100;

const int delta_t = (t_max - t_min)/(double)n_steps;

double last_percent(-1.);

class THit{
public:
	int col;
	int row;
	double x;
	double y;
	double z;
	int time;
	int event;
	int plane;
  int MCtrk;
  THit():col(0),row(0), x(0), y(0), z(0), time(0), event(0), plane(0), MCtrk(-1){
		;
	}
};

typedef vector<THit> THits;

THits generate_hits();

class TCell{ // stores index for a hit upstream and downstream
public:
	int ihit_upstream;
	int ihit_downstream;
	double dxdz; // direction of the cell in xz plane
	double dydz; // direction of the cell in yz plane
};

typedef vector<TCell> TCells;

class TBreakingAngle{
public:
	int icell_upstream;
	int icell_downstream;
	double angle_dxdz;
	double angle_dydz;
};

typedef vector<TBreakingAngle> TBreakingAngles;

class TTrack{
public:
	//THits hits; // hits assigned to the track
	TCells cells; // cells indexing the hits above
	TBreakingAngles breaking_angles; // breaking angles between two cells;
};

typedef vector<TTrack> TTracks;

// search track candidates in a simple cellular automaton style
// assuming full efficiency of all planes
TTracks SearchTracks(const THits& hits){
	TTracks result;

	// build cells between all possible combinations of hits in different planes
	TCells cells;
	for (int ihit0 = 0; ihit0 < hits.size(); ihit0++){
		for (int ihit1 = ihit0+1; ihit1 < hits.size(); ihit1++){
			//if (ihit0 == ihit1) continue;
			int hit_up = ihit0;
			int hit_do = ihit1;
			// use the correct order
			if (hits[hit_up].plane > hits[hit_do].plane){
				hit_up = ihit1;
				hit_do = ihit0;
			}
			// cut neighboring planes only
			// (no missing planes => != 1)
			// missing planes => <= 0
			// in both cases:
			// no hits in the same plane
			if ((hits[hit_do].plane - hits[hit_up].plane) <= 0){
				continue;
			}
			//cout << hits[hit_do].plane << " " << hits[hit_up].plane << endl;
			// build cells
			TCell cell;
			cell.ihit_downstream = hit_do;
			cell.ihit_upstream = hit_up;
			double dz = (hits[hit_do].z - hits[hit_up].z);
			cell.dxdz = (hits[hit_do].x - hits[hit_up].x) / dz;
			cell.dydz = (hits[hit_do].y - hits[hit_up].y) / dz;
			//cout << " angles are " << dz << " " << cell.dxdz << " " << cell.dydz << endl;
			cells.push_back(cell);
			//cout << hits[cells[cells.size()-1].ihit_downstream].plane << " " << hits[cells[cells.size()-1].ihit_upstream].plane << endl << endl;
		}
	}
	// connect cells to each other and calculate the breaking angles
	// use only neighboring cells
	// I think this algorithm will currently produce all cell connections twice
	// to be rechecked!!!
	TBreakingAngles breaking_angles;
	TCells trackcells;
	//cout << " running over " << cells.size() << " cells " << endl;
	for (int icell0 = 0; icell0 < cells.size(); icell0++){
		for (int icell1 = icell0+1; icell1 < cells.size(); icell1++){
			// skip same cells
			//if (icell0 == icell1) continue;
			//cout << " searching cell combo " << icell1 << " " << icell0 << endl;
			int icellup = 0;
			int icelldo = 0;
			// use the correct order and
			// connect only cells sharing hits
			//cout << hits[cells[icell0].ihit_downstream].plane << " " << hits[cells[icell1].ihit_upstream].plane << endl;
			if (cells[icell0].ihit_downstream == cells[icell1].ihit_upstream){
				icellup = icell0;
				icelldo = icell1;
			} else {
				if (cells[icell0].ihit_upstream == cells[icell1].ihit_downstream){
					icellup = icell1;
					icelldo = icell0;
				} else {
					continue;
				}
			}
			//cout << cells[icell0].ihit_downstream << " " << cells[icell1].ihit_upstream << endl;
			// connect cells and calculate the breaking angles
			TBreakingAngle breaking_angle;
			breaking_angle.icell_upstream = trackcells.size();
			trackcells.push_back(cells[icellup]);
			breaking_angle.icell_downstream = trackcells.size();
			trackcells.push_back(cells[icelldo]);
			double dxdz_up = cells[icellup].dxdz;
			double dydz_up = cells[icellup].dydz;
			double dxdz_do = cells[icelldo].dxdz;
			double dydz_do = cells[icelldo].dydz;
			// calculate the angle difference, keep the sign
			//if (dxdz_up > dxdz_do)
				breaking_angle.angle_dxdz = dxdz_up - dxdz_do;
			//else
			//	breaking_angle.angle_dxdz = dxdz_do - dxdz_up;
			//if (dydz_up > dydz_do)
				breaking_angle.angle_dydz = dydz_up - dydz_do;
			//else
			//	breaking_angle.angle_dydz = dydz_do - dydz_up;
			//cout << " breaking angles are " << breaking_angle.angle_dxdz << " and " << breaking_angle.angle_dydz << endl;
			breaking_angles.push_back(breaking_angle);
			// create a track segment
			TTrack track;
			track.breaking_angles = breaking_angles;
			track.cells = trackcells;
			//track.hits =
			result.push_back(track);
			breaking_angles.clear();
			trackcells.clear();
		}
	}

	return result;
}

void DrawProgressBar(int len, double percent) {
	if ((int) (last_percent * 100) == (int) (percent * 100))
		return;
	//cout << " drawing " << endl;
	cout << "\x1B[2K"; // Erase the entire current line.
	cout << "\x1B[0E"; // Move to the beginning of the current line.
	string progress;
	for (int i = 0; i < len; ++i) {
		if (i < static_cast<int> (len * percent)) {
			progress += "=";
		} else {
			progress += " ";
		}
	}
	cout << "[" << progress << "] " << (static_cast<int> (100 * percent))
			<< "%";
	flush( cout); // Required.
	last_percent = percent;
}

void track_finder(){
	TCanvas canvas_distributions("canvas_distributions", "distributions", 800, 800);
	TH2F* hist_XY[4];
	hist_XY[0] = new TH2F("hist_XY_plane_0", "xy distribution at plane 0", 49, -2e-3, 2e-3, 49, -2e-3, 2e-3);
	hist_XY[1] = new TH2F("hist_XY_plane_1", "xy distribution at plane 1", 49, -2e-3, 2e-3, 49, -2e-3, 2e-3);
	hist_XY[2] = new TH2F("hist_XY_plane_2", "xy distribution at plane 2", 49, -2e-3, 2e-3, 49, -2e-3, 2e-3);
	hist_XY[3] = new TH2F("hist_XY_plane_3", "xy distribution at plane 3", 49, -2e-3, 2e-3, 49, -2e-3, 2e-3);

	TH2F* hist_dxdz[2]; // breaking angles between planes
	hist_dxdz[0] = new TH2F("hist_dxdz_plane_012", "dx/dz and dy/dz distribution between plane 0,1 and 2", 49, -2e-2, 2e-2, 49, -2e-2, 2e-2);
	hist_dxdz[1] = new TH2F("hist_dxdz_plane_123", "dx/dz and dy/dz distribution between plane 1,2 and 3", 49, -2e-2, 2e-2, 49, -2e-2, 2e-2);
//	hist_dxdz[2] = new TH2F("hist_dxdz_plane_23", "dxdz distribution between plane 2 and 3", 49, -2e-2, 2e-2, 49, -2e-2, 2e-2);

	TH2F* hist_dxdz_track; // breaking angles between planes
	hist_dxdz_track = new TH2F("hist_dxdz_track", "dx/dz and dy/dz direction distribution", 49, -0.02, 0.02, 49, -0.02, 0.02);

	TH2I *hist_trk_mult = new TH2I("hist_trk_mult",";MC trks/ev; REC trks/ev",20,0,20,20,0,20);
	int mean_hits(0);
	int mean_tracks(0);

	int nevents = 1000000;
	for (int ievent = 0; ievent < nevents; ievent++){
		DrawProgressBar(50, (ievent + 1) / ((double) nevents));
		THits hits = generate_hits();
		for (auto ihit = hits.begin(); ihit != hits.end(); ihit++){
			hist_XY[ihit->plane]->Fill(ihit->x, ihit->y);
		}
		mean_hits += hits.size();
		int MCMCtrks=0;//number of MC trks left hits in LMD
		if((hits.size())>1){
		  cout<<"------ Event "<< ievent<<" ------- "<<endl;
		  cout<<"MC MC ids: ";
		  int MCMCid=-1;
		  for (auto ihit = hits.begin(); ihit != hits.end(); ihit++){
		    int MCtrk_cur = ihit->MCtrk;
		    if(MCMCid!=MCtrk_cur){
		      MCMCid=MCtrk_cur;
		      MCMCtrks++;
		  }

		  cout<<" "<<MCtrk_cur;
		}
		cout<<endl;
		}

		TTracks tracks = SearchTracks(hits);
		mean_tracks += tracks.size();
		if(tracks.size()>0) hist_trk_mult->Fill(MCMCtrks,tracks.size());
		for (auto itrack = tracks.begin(); itrack != tracks.end(); itrack++){
			double dxdz_mean(0);
			double dydz_mean(0);
			for (auto icell = itrack->cells.begin(); icell != itrack->cells.end(); icell++){
				dxdz_mean += (*icell).dxdz;
				dydz_mean += (*icell).dydz;
			}
			dxdz_mean /= itrack->cells.size();
			dydz_mean /= itrack->cells.size();
			hist_dxdz_track->Fill(dxdz_mean, dydz_mean);
			for (auto iangle = itrack->breaking_angles.begin(); iangle != itrack->breaking_angles.end(); iangle++){
				// find the starting plane of connected cells
				int iplane_start = (hits[itrack->cells[(iangle->icell_upstream)].ihit_upstream]).plane;
				//cout << " breaking angles are " << iangle->angle_dxdz << " and " << iangle->angle_dydz << endl;
				hist_dxdz[iplane_start]->Fill(iangle->angle_dxdz, iangle->angle_dydz);
				//hist_dydz[iplane_start]->Fill(iangle->angle_dydz);
			}
		}
	}

  	cout << endl << " some statistics for " << nevents << " events " << endl;
	cout << " total number of hits processed: " << mean_hits << endl;
	cout << " total number of tracks found: " << mean_tracks << endl;
	cout << " ratio tracks / hits " << ((double) mean_tracks)/((double) mean_hits) << endl;
	cout << " mean number of hits " << ((double) mean_hits)/((double) nevents) << endl;
	cout << " mean number of tracks " << ((double) mean_tracks)/((double) nevents) << endl;

	canvas_distributions.Divide(2,2);
	for (int iplane = 0; iplane < 4; iplane++){
		canvas_distributions.cd(iplane+1);
		hist_XY[iplane]->Draw("colz");
	}
	canvas_distributions.Print("hit_ditros.pdf(");
	canvas_distributions.cd(3);
	hist_dxdz_track->Draw("colz");
	canvas_distributions.cd(4);
	gPad->Clear();
        //hist_trk_mult->Draw("colz");
	//canvas_distributions.Print("hit_ditros.pdf(");
	for (int iplane = 0; iplane < 2; iplane++){
		canvas_distributions.cd(iplane+1);
		hist_dxdz[iplane]->Draw("colz");
	}
	canvas_distributions.Print("hit_ditros.pdf(");
	canvas_distributions.cd(1);
	hist_trk_mult->Draw("colz");
	canvas_distributions.cd(2);
	gPad->Clear();
	canvas_distributions.cd(3);
	gPad->Clear();
	canvas_distributions.cd(4);
	gPad->Clear();
	canvas_distributions.Print("hit_ditros.pdf)");
}

double stepfunction(const double* x, const double* pixelsize){
	/*double result = x[0];
	result *= 1e6;
	result /= 80;
	result = floor(result);
	result *= 80.;
	result *= 1e-6;
	return result;*/
	return (floor((x[0]*1e6+pixelsize[0]/2)/pixelsize[0])*pixelsize[0])*1e-6;
}

// simulate electrons going through 4 HV-MAPS mupix6 prototypes
// behind the tagger (setup of march 2015 at MAMI)
// simulate ntracks which hits are distributed over nevents
THits generate_hits(){
	const int ntracks_mean = 5;
	int ntracks = gRandom->Poisson(ntracks_mean);
	const int nevents = 1; // not used yet
	const double momentum = 0.9*1.5; // GeV/c
	// direction from the focal plane. 0 is perpendicular to the detector plane
	// unit is rad
	const double dxdz_mean = -0.005;//-0.01;//-0.1;
	const double dydz_mean = +0.003;
	const double dxdz_sigma = 0.001;
	const double dydz_sigma = 0.001;
	const double posx = 0.1; // location of the starting point +/- [m]
	const double posy = 0.01; // location of the starting point +/- [m]
	THits result;
	const double pixelsize = 80; // um
	// z positions of planes in m
	double z[5] = {0.0, 0.9, 1.0, 1.1, 1.2};
		// acceptance in the plane
	const double xmin = -1.5e-3;
	const double xmax = 1.5e-3;

	const double beta = 1/sqrt(1+pow(0.511/momentum,2));
	//cout << "the beta factor for " << momentum << " GeV/c electrons is " << beta << endl;
	const double x_over_x0 = 2.329/21.82*350/1000/10;
	//cout << "x/X0 is" << x_over_x0 << endl;
	const double theta_0_slicon = 13.6/beta/(momentum*1e3)*sqrt(x_over_x0)*(1+0.038*log(x_over_x0));
	//cout << "scattering angle of 350 mum silicon is " << theta_0_slicon << endl;

	double dxdz =  gRandom->Gaus(dxdz_mean, dxdz_sigma);
	double dydz =  gRandom->Gaus(dydz_mean, dydz_sigma);

	TVector3 dir(0., 0., momentum);
	dir.SetXYZ(dxdz, dydz, momentum);

	TVector3 pos(gRandom->Uniform(-posx, posx), gRandom->Uniform(-posy, posy), 0.); // m

	// propagate to the first plane
	for (int itrack = 0; itrack < ntracks; itrack++)
	for (int iplane = 0; iplane < 4; iplane++){
		pos += dir.Unit()*(z[iplane+1]-z[iplane]);
		TVector3 pos_plane_0 = pos;
		double X = pos.X();
		double Y = pos.Y();
		// simulate a granularity of the detector in theta of XX mum
		//		if (xmin < X && X < xmax && xmin < Y && Y < xmax)
		{
			TVector3 pos_reco_plane(stepfunction(&X,&pixelsize),
					stepfunction(&Y,&pixelsize),
					z[iplane+1]);
			THit hit;
			hit.x = pos_reco_plane.X();
			hit.y = pos_reco_plane.Y();
			hit.z = pos_reco_plane.Z();
			hit.plane = iplane;
			hit.MCtrk = itrack;
			// missalign
			if (iplane == 1){
				double _x = hit.x;
				double _y = hit.y;
				hit.x = _x *  cos(0.06) + _y * sin(0.06) + 0.00016;
				hit.y = _x * -sin(0.06) + _y * cos(0.06) - 0.00004;
			}

			if (iplane == 2){
				double _x = hit.x;
				double _y = hit.y;
				hit.x = _x *  cos(0.1) + _y * -sin(0.1) - 0.00020;
				hit.y = _x *  sin(0.1) + _y *  cos(0.1) - 0.00012;
			}

			if (iplane == 3){
				double _x = hit.x;
				double _y = hit.y;
				hit.x = _x *  cos(0.04) + _y * -sin(0.04) + 0.00012;
				hit.y = _x *  sin(0.04) + _y *  cos(0.04) + 0.00020;
			}

			//if (iplane != 3)
			if (xmin < X && X < xmax && xmin < Y && Y < xmax){
				result.push_back(hit);
			}

			// smear angle due to silicon in the xz plane
			double tan_theta_new = tan(atan(dir.X()/dir.Z())+gRandom->Gaus(0,theta_0_slicon));
			double l_x = dir.X()*dir.X()+dir.Z()*dir.Z();
			double x_new = sqrt(l_x / (1/(tan_theta_new*tan_theta_new)+1));
			// recover the lost sign
			if (tan_theta_new < 0 && x_new > 0) x_new = -x_new;
			double z_new = sqrt(l_x - x_new*x_new);
			dir.SetX(x_new);
			dir.SetZ(z_new);
			// same smearing procedure in the yz plane
			tan_theta_new = tan(atan(dir.Y()/dir.Z())+gRandom->Gaus(0,theta_0_slicon));
			double l_y = dir.Y()*dir.Y()+dir.Z()*dir.Z();
			double y_new = sqrt(l_y / (1/(tan_theta_new*tan_theta_new)+1));
			// recover the lost sign
			if (tan_theta_new < 0 && y_new > 0) y_new = -y_new;
			z_new = sqrt(l_y - y_new*y_new);
			dir.SetY(y_new);
			dir.SetZ(z_new);
		} //else {
		//	break;
		//}
	}
	return result;
}


void t_to_theta_transformation(){
	for (int istep = 0; istep < n_steps; istep++){
		double t = t_min + istep*delta_t;
		TLorentzVector p_in;
		TLorentzVector p_out;
		TLorentzVector p_bar_in;
		TLorentzVector p_bar_out;

		double momentum_in = sqrt(energy*energy - mass*mass);
		p_bar_in.SetPxPyPzE(0., 0., momentum_in, energy); // incoming beam
		p_in.    SetPxPyPzE(0., 0., 0., mass); // target proton at rest
	}
}

vector<unsigned int> evttimes; // in ns
vector<unsigned int> evtnbrs;
vector<double> evttimes_graph; // in ns
vector<double> evtnbrs_graph;
map<unsigned int,unsigned int> evtframes; // time on a basis of 25ns with a counter for the number of events in this frame

void timebased_study(){
	TRandom3* rnd = new TRandom3();
	double cross_section = 15.67e-27 * 2.; // [cm^-2]
	double luminosity = 2e32; // [cm^-2 s^-1]
	unsigned int nEvents = cross_section*luminosity; // per second
	double timelength = 1e9; // in [ns] 0.8 s for a reference luminosity of 1e32 cm^-2 s^-1
	cout << " simulating " << nEvents << " over a period of " << timelength*1e-9 << " s" << endl;
	for (unsigned int iEvent = 0; iEvent < nEvents; iEvent++){
		DrawProgressBar(50, (iEvent + 1) / ((double) nEvents));
		double time = rnd->Uniform(timelength);
		unsigned int time1ns = (unsigned int) floor(time);
		unsigned int time25ns = time/25;
		evttimes.push_back(time1ns);
		evtnbrs.push_back(iEvent);
		evttimes_graph.push_back(time);
		evtnbrs_graph.push_back((double) iEvent);
		evtframes[time25ns]++;
	}
	cout << endl << " analyzing multiplicities " << endl;
	TH1D* hist_evtframe_multiplicity = new TH1D("hist_evtframe_multiplicity", "multiplicity in event frames", 100, -0.5, 99.5);
	// go through all time frames in a period of <timelength> s
	unsigned int nFrames = (unsigned int) floor(timelength/25.);
	for (unsigned int iFrame = 0; iFrame < nFrames; iFrame++){
		DrawProgressBar(50, (iFrame + 1) / ((double) nFrames));
		hist_evtframe_multiplicity->Fill(evtframes[iFrame]);
	}
	TCanvas* canvas_multiplicity = new TCanvas("canvas_multiplicity", "multiplicity distributions", 600, 600);
	hist_evtframe_multiplicity->Draw("hist text");
	gPad->Update();

	cout << endl << " viewing histogram numbers " << endl;

	unsigned int entries0 = hist_evtframe_multiplicity->GetBinContent(1);
	unsigned int entries1 = hist_evtframe_multiplicity->GetBinContent(2);
	unsigned int entries_meanw0(0);
	unsigned int entries_meanwo0(0);
	unsigned int entries_more1;
	unsigned int frames_all(0);
	unsigned int frames_allwo0(0);
	for (int ibin = 1; ibin < 100; ibin++){
		entries_meanw0 += hist_evtframe_multiplicity->GetBinContent(ibin)*(ibin-1);
		frames_all += hist_evtframe_multiplicity->GetBinContent(ibin);
		if (ibin > 1){
			entries_meanwo0 += hist_evtframe_multiplicity->GetBinContent(ibin)*(ibin-1);
			frames_allwo0 += hist_evtframe_multiplicity->GetBinContent(ibin);
		}
		if (ibin > 2)
			entries_more1 += hist_evtframe_multiplicity->GetBinContent(ibin);
	}
	cout << frames_all << " time frames with a 25 ns basis found in the histogram " << endl;
	cout << " we should expect " << nFrames << " time frames " << endl;
	cout << entries0 << " time frames with no hits or " << (double)entries0/(double)frames_all*100. << "% of all frames" << endl;
	cout << entries1 << " time frames with one hit or " << (double)entries1/(double)frames_all*100. << "% of all frames" << endl;
	cout << " corresponding to " << (double)entries1/(double)frames_allwo0*100. << "% of all frames without 0 ones" << endl;
	cout << entries_more1 << " time frames with more than 1 hit or " << (double)entries_more1/(double)frames_all*100. << "% of all events" << endl;
	cout << " corresponding to " << (double)entries_more1/(double)frames_allwo0*100. << "% of all events without 0 entries" << endl;
	cout << " mean number of hits per frame is " << (double) entries_meanw0 / (double) frames_all << endl;
	cout << " mean number of hits per frame without 0 is " << (double) entries_meanwo0 / (double) frames_allwo0 << endl;

	cout << endl << " sorting " << endl;
	sort(evttimes_graph.begin(), evttimes_graph.end());
	cout << endl << " calculating some numbers " << endl;
	double mean_time(0);
	double min_time = evttimes_graph[0];
	nEvents = evttimes_graph.size();
	double max_time = evttimes_graph[nEvents-1];

	cout << endl << " simulated " << nEvents << " events in a time window between " << min_time*1e-9 << " s and " << max_time*1e-9 << " s " << endl;

	for (unsigned int iEvent = 1; iEvent < nEvents; iEvent++){
		DrawProgressBar(50, (iEvent + 1) / ((double) nEvents));
		double delta_t = evttimes_graph[iEvent]-evttimes_graph[iEvent-1];
		if (delta_t < 0) cout << " warning: vector is not sorted! " << endl;
		mean_time+=delta_t;
	}
	mean_time /= (double) nEvents;
	cout << endl << " mean time between 2 consecutive events is " << mean_time << " ns " << endl;
	cout << " what corresponds to a mean rate of " << 1/mean_time*1e3 << " MHz " << endl;
	cout << " this should fit the mean number of " << nEvents / timelength*1e9 << " events / s " << endl;
	//cout << " drawing " << endl;
	//TGraph* graph = new TGraph(nEvents, &evttimes_graph[0], &evtnbrs_graph[0]);
	//graph->Draw("ALP");

	/* results are
	 *
	 * ./Debug/sandbox
 simulating 6268000 over a period of 1 s
[==================================================] 100%
 analyzing multiplicities
[==================================================] 100%
 viewing histogram numbers
40000000 time frames with a 25 ns basis found in the histogram
 we should expect 40000000 time frames
34197585 time frames with no hits or 85.494% of all frames
5360628 time frames with one hit or 13.4016% of all frames
 corresponding to 92.3862% of all frames without 0 ones
474554 time frames with more than 1 hit or 1.18639% of all events
 corresponding to 8.17856% of all events without 0 entries
 mean number of hits per frame is 0.1567
 mean number of hits per frame without 0 is 1.08024

 sorting

 calculating some numbers

 simulated 6268000 events in a time window between 4.16068e-07 s and 0.999999 s
[==================================================] 100%
 mean time between 2 consecutive events is 159.54 ns
 what corresponds to a mean rate of 6.26801 MHz
 this should fit the mean number of 6.268e+06 events / s
	 *
	 */
}

void MickeyMouse_MC_asymetries(){
	const int nmomenta = 1;//3;
	double momenta[nmomenta] = {15};//{1.5, 4.06, 15.00};
	double pixelsize = 80; // um
	double ip_size = 1.; // mm
	double ip_div = 1.; // mrad
	double ip_width_z = 5.; // mm
	// z positions of planes in m
	vector < double > z; z.push_back(11.5); z.push_back(11.7); z.push_back(11.8); z.push_back(11.9);
	// acceptance in the plane (already the power of 2)
	const double Rmin = 0.035*0.035;
	const double Rmax = (0.035+0.06)*(0.035+0.06);
	// cut for a cone starting from ip
	// angles for a cone cut
	const double theta_cut_low = 4.; // mrad
	const double theta_cut_high = 7.; // mrad
	vector < double > Rmin_cuts;
	vector < double > Rmax_cuts;
	for (int iplane = 0; iplane < 4; iplane++ ){
		Rmin_cuts.push_back(theta_cut_low*1e-3*z[iplane]*theta_cut_low*1e-3*z[iplane]);
		Rmax_cuts.push_back(theta_cut_high*1e-3*z[iplane]*theta_cut_high*1e-3*z[iplane]);
	}

	TCanvas canvas("canvas", "canvas", 600, 1000);
	canvas.Divide(2,1);
	canvas.cd(1);
	TCanvas canvas_distributions("canvas_distributions", "distributions", 800, 800);
	canvas_distributions.Divide(2,2);

	//TFile treefile("asymmetries_results.root", "RECREATE");
	TTree results("asymmetry_results", "Asymmetry results");
	double offsetx, offsety, tiltx, tilty, momentum,
		offsetx_fit, offsety_fit, tiltx_fit, tilty_fit,
		meanx_0, meany_0, rmsx_0, rmsy_0,
		meanx_1, meany_1, rmsx_1, rmsy_1,
		meanx_2, meany_2, rmsx_2, rmsy_2,
		meanx_3, meany_3, rmsx_3, rmsy_3;
	int	count_0[10] = {0,0,0,0,0,0,0,0,0,0};
	int count_1[10] = {0,0,0,0,0,0,0,0,0,0};
	int count_2[10] = {0,0,0,0,0,0,0,0,0,0};
	int count_3[10] = {0,0,0,0,0,0,0,0,0,0};
	TLine line;

	results.Branch("offsetx", &offsetx);
	results.Branch("offsety", &offsety);
	results.Branch("tiltx", &tiltx);
	results.Branch("tilty", &tilty);
	results.Branch("momentum", &momentum);
	results.Branch("offsetx_fit", &offsetx_fit);
	results.Branch("offsety_fit", &offsety_fit);
	results.Branch("tiltx_fit", &tiltx_fit);
	results.Branch("tilty_fit", &tilty_fit);
	results.Branch("meanx_0", &meanx_0);
	results.Branch("meany_0", &meany_0);
	results.Branch("rmsx_0", &rmsx_0);
	results.Branch("rmsy_0", &rmsy_0);
	results.Branch("meanx_1", &meanx_1);
	results.Branch("meany_1", &meany_1);
	results.Branch("rmsx_1", &rmsx_1);
	results.Branch("rmsy_1", &rmsy_1);
	results.Branch("meanx_2", &meanx_2);
	results.Branch("meany_2", &meany_2);
	results.Branch("rmsx_2", &rmsx_2);
	results.Branch("rmsy_2", &rmsy_2);
	results.Branch("meanx_3", &meanx_3);
	results.Branch("meany_3", &meany_3);
	results.Branch("rmsx_3", &rmsx_3);
	results.Branch("rmsy_3", &rmsy_3);
	results.Branch("count_0", count_0, "count_0[10]/I");
	results.Branch("count_1", count_1, "count_1[10]/I");
	results.Branch("count_2", count_2, "count_2[10]/I");
	results.Branch("count_3", count_3, "count_3[10]/I");

	// calculate the number of loops
	int nloops = 0;
	for (offsetx = -2.; offsetx <= 2.; offsetx+=2.) // mm
	for (offsety = -2.; offsety <= 2.; offsety+=2.)
	for (tiltx = -.5; tiltx <= .5; tiltx+=.5) // mrad
	for (tilty = -.5; tilty <= .5; tilty+=.5) // mrad
	for (int imom = 0; imom < nmomenta; imom++){
		nloops++;
	}

	int iloop = 0;

	TGraphAsymmErrors* graph_model = NULL;
	TH1* model_distribution = NULL;
	for (int imom = 0; imom < nmomenta; imom++){
		canvas.cd(1);
		momentum = momenta[imom];
		// try to retrieve the DPM density function to model the angular distribution
		stringstream modelfilename;
		modelfilename << "DPMModels_" << momentum << "_model.root";
		TFile model(modelfilename.str().c_str(), "READ");
		delete graph_model;
		TGraphAsymmErrors* graph_model = (TGraphAsymmErrors*) model.Get("Graph");
		TH1F model_distribution("model_distribution", "distribution", 1000, 1., 10.);
		model.Close();
		if (!graph_model){
			cout << " Warning: no density function found in " << modelfilename.str() << endl;
		} else {
			//graph_model->Draw("ALP");
			//gPad->Update();
			for (double ibin = 1.; ibin <= model_distribution.GetXaxis()->GetNbins(); ibin++){
				double theta = model_distribution.GetXaxis()->GetBinCenter(ibin); // mrad
				model_distribution.Fill(theta, graph_model->Eval(theta));
			}
			gPad->Clear();
			model_distribution.Draw();
			gPad->Update();
		}
	for (offsetx = -2.; offsetx <= 2.; offsetx+=2.){ // mm
	for (offsety = -2.; offsety <= 2.; offsety+=2.){
	for (tiltx = -.5; tiltx <= .5; tiltx+=.5){ // mrad
	for (tilty = -.5; tilty <= .5; tilty+=.5){ // mrad
		iloop++;
		cout << " **** " << iloop << " of " << nloops << " ******** " << endl;
		cout << " **** " << iloop*100/nloops << "% **********" << endl;
		stringstream filename;
		filename << "asymmetry_studies_"<< momentum << "_GeVc_x" << offsetx << "_mm_y_" << offsety << "_mm_dxdz_" << tiltx << "_mrad_dydz_" << tilty << "_mrad";
		//TFile file((filename.str()+".root").c_str(), "RECREATE");

		TH2F hist_gen_theta_phi ("hist_gen_theta_phi", "theta phi generated distribution", 100, 2e-3, 9e-3, 100, -3.141, 3.141);
		TH2F hist_gen_x_y ("hist_gen_x_y", "x y generated distribution", 100, -5e-3, 5e-3, 100, -5e-3, 5e-3);
		TH2F hist_gen_dx_dy ("hist_gen_dx_dy", "dx/dz dy/dz generated distribution", 100, -15e-3, 15e-3, 100, -15e-3, 15e-3);
		TH2F hist_XY_plane_0 ("hist_XY_plane_0", "xy distribution at plane 0", 100, -0.100, 0.100, 100, -0.100, 0.100);
		TH2F hist_XY_plane_1 ("hist_XY_plane_1", "xy distribution at plane 1", 100, -0.100, 0.100, 100, -0.100, 0.100);
		TH2F hist_XY_plane_2 ("hist_XY_plane_2", "xy distribution at plane 2", 100, -0.100, 0.100, 100, -0.100, 0.100);
		TH2F hist_XY_plane_3 ("hist_XY_plane_3", "xy distribution at plane 3", 100, -0.100, 0.100, 100, -0.100, 0.100);
		//TH2F hist_phi_count("hist_phi_count", "counts over plane number", 4, -0.5, 3.5, 10, -3.141, 3.141);
		//TH2F hist_center_of_gravity("hist_center_of_gravity", "center of gravity over plane number", 4, -0.5, 3.5, 10, -3.141, 3.141);

		// calculate some multiple scattering constants

		double beta = 1/sqrt(1+pow(0.938/momentum,2));
		cout << "the beta factor for " << momentum << " GeV/c antiprotons is " << beta << endl;
		double x_over_x0 = 2.329/21.82*350/1000/10;
		cout << "x/X0 is" << x_over_x0 << endl;
		double theta_0_slicon = 13.6/beta/(momentum*1e3)*sqrt(x_over_x0)*(1+0.038*log(x_over_x0));
		cout << "scattering angle of 350 mum silicon is " << theta_0_slicon << endl;
		// assuming a multiple scattering effect in the transition region of a silicon plane of half thickness
		double theta_0_kapton = 13.6/beta/(momentum*1e3)*sqrt(x_over_x0/2.)*(1+0.038*log(x_over_x0/2.));
		for (int ievent = 0; ievent < 1000000*1.7; ievent++){
			double phi =  gRandom->Uniform(-3.141, 3.141);
			double theta;
			if (graph_model) theta = model_distribution.GetRandom()*1e-3;
			else theta = gRandom->Uniform(1e-3, 10e-3);
			TVector3 dir(0., 0., momentum);
			dir.SetTheta(theta);
			dir.SetPhi(phi);
			TVector3 pos(gRandom->Gaus(offsetx, ip_size)*1e-3,gRandom->Gaus(offsety, ip_size)*1e-3,gRandom->Uniform(-ip_width_z, ip_width_z)*1e-3); // m
			// rotate according to a tilted and smeared beam
			TVector3 dir_beam(gRandom->Gaus(tiltx, ip_div)*1e-3,gRandom->Gaus(tilty, ip_div)*1e-3,1.);
			//cout << "dir_beam " << dir_beam.X()/dir_beam.Z()-tiltx*1e-3 << " " << dir_beam.Y()/dir_beam.Z()-tilty*1e-3 << endl;
			dir.RotateUz(dir_beam.Unit());
			//cout << "dir " << (dir.X()/dir.Z()*dir.X()/dir.Z()+dir.Y()/dir.Z()*dir.Y()/dir.Z()) tiltx*1e-3 << " " << dir.Y()/dir.Z()-tilty*1e-3 << endl;
			hist_gen_x_y.Fill(pos.X(), pos.Y());
			hist_gen_theta_phi.Fill(dir.Theta(), dir.Phi());
			hist_gen_dx_dy.Fill(dir.X()/dir.Z(), dir.Y()/dir.Z());

			//hist_gen.Fill(dir.Theta(), dir.Phi());
			// propagate to the transition region at about 11m
			pos += dir.Unit()*11.; // sure propagation is not fully correct as it is not propagated to a plane
			// smear angle due to silicon in the xz plane
			double tan_theta_new = tan(atan(dir.X()/dir.Z())+gRandom->Gaus(0,theta_0_kapton));
			double l_x2 = dir.X()*dir.X()+dir.Z()*dir.Z();
			double x_new = sqrt(l_x2 / (1/(tan_theta_new*tan_theta_new)+1));
			double z_new = sqrt(l_x2 - x_new*x_new);
			// recover the lost sign
			if (dir.X() < 0 && x_new > 0) x_new = -x_new;
			dir.SetX(x_new);
			dir.SetZ(z_new);
			// same smearing procedure in the yz plane
			tan_theta_new = tan(atan(dir.Y()/dir.Z())+gRandom->Gaus(0,theta_0_kapton));
			double l_y2 = dir.Y()*dir.Y()+dir.Z()*dir.Z();
			double y_new = sqrt(l_y2 / (1/(tan_theta_new*tan_theta_new)+1));
			// recover the lost sign
			if (dir.Y() < 0 && y_new > 0) y_new = -y_new;
			z_new = sqrt(l_y2 - y_new*y_new);
			dir.SetY(y_new);
			dir.SetZ(z_new);
			// propagate to the first plane
			pos += dir.Unit()*0.5;
			TVector3 pos_plane_0 = pos;
			//hist_phi_count.Fill(0, pos_plane_0.Phi());
			// simulate a granularity of the detector in theta of XX mum
			double X = pos.X();
			double Y = pos.Y();
			double RR = X*X+Y*Y;
			if (Rmin_cuts[0] < RR && RR < Rmax_cuts[0]){//(Rmin < RR && RR < Rmax){
				TVector3 pos_reco_plane_0(stepfunction(&X,&pixelsize),
						stepfunction(&Y,&pixelsize),
						11.500);
				int imodule = (int) floor(((pos_reco_plane_0.Phi()+3.1415)/(2*3.1415))*10.);
				if (imodule > 9){ /*cout << imodule << endl;*/ imodule = 0; }
				if (imodule < 0){ cout << imodule << endl; imodule = 0; }
				count_0[imodule]++;
				if (imodule == 0 || imodule == 3 || imodule == 5 || imodule == 8)
				hist_XY_plane_0.Fill(pos_reco_plane_0.X(), pos_reco_plane_0.Y());
			}
			// ****** plane 1 ********
			// smear angle due to silicon in the xz plane
			tan_theta_new = tan(atan(dir.X()/dir.Z())+gRandom->Gaus(0,theta_0_slicon));
			l_x2 = dir.X()*dir.X()+dir.Z()*dir.Z();
			x_new = sqrt(l_x2 / (1/(tan_theta_new*tan_theta_new)+1));
			// recover the lost sign
			if (dir.X() < 0 && x_new > 0) x_new = -x_new;
			z_new = sqrt(l_x2 - x_new*x_new);
			dir.SetX(x_new);
			dir.SetZ(z_new);
			// same smearing procedure in the yz plane
			tan_theta_new = tan(atan(dir.Y()/dir.Z())+gRandom->Gaus(0,theta_0_slicon));
			l_y2 = dir.Y()*dir.Y()+dir.Z()*dir.Z();
			y_new = sqrt(l_y2 / (1/(tan_theta_new*tan_theta_new)+1));
			// recover the lost sign
			if (dir.Y() < 0 && y_new > 0) y_new = -y_new;
			z_new = sqrt(l_y2 - y_new*y_new);
			dir.SetY(y_new);
			dir.SetZ(z_new);
			// propagate to the second plane
			pos += dir.Unit()*0.2;
			TVector3 pos_plane_1 = pos;
			X = pos.X();
			Y = pos.Y();
			RR = X*X+Y*Y;
			if (Rmin_cuts[1] < RR && RR < Rmax_cuts[1]){//(Rmin < RR && RR < Rmax){
				TVector3 pos_reco_plane_1(stepfunction(&X,&pixelsize),
						stepfunction(&Y,&pixelsize),
						11.700);
				int imodule = (int) floor(((pos_reco_plane_1.Phi()+3.1415)/(2*3.1415))*10.);
				if (imodule > 9){ /*cout << imodule << endl;*/ imodule = 0; }
				if (imodule < 0){ cout << imodule << endl; imodule = 0; }
				count_1[imodule]++;
				if (imodule == 0 || imodule == 3 || imodule == 5 || imodule == 8)
					hist_XY_plane_1.Fill(pos_reco_plane_1.X(), pos_reco_plane_1.Y());
			}
			// ****** plane 2 ********
			// smear angle due to silicon in the xz plane
			tan_theta_new = tan(atan(dir.X()/dir.Z())+gRandom->Gaus(0,theta_0_slicon));
			l_x2 = dir.X()*dir.X()+dir.Z()*dir.Z();
			x_new = sqrt(l_x2 / (1/(tan_theta_new*tan_theta_new)+1));
			// recover the lost sign
			if (dir.X() < 0 && x_new > 0) x_new = -x_new;
			z_new = sqrt(l_x2 - x_new*x_new);
			dir.SetX(x_new);
			dir.SetZ(z_new);
			// same smearing procedure in the yz plane
			tan_theta_new = tan(atan(dir.Y()/dir.Z())+gRandom->Gaus(0,theta_0_slicon));
			l_y2 = dir.Y()*dir.Y()+dir.Z()*dir.Z();
			y_new = sqrt(l_y2 / (1/(tan_theta_new*tan_theta_new)+1));
			// recover the lost sign
			if (dir.Y() < 0 && y_new > 0) y_new = -y_new;
			z_new = sqrt(l_y2 - y_new*y_new);
			dir.SetY(y_new);
			dir.SetZ(z_new);
			// propagate to the second plane
			pos += dir.Unit()*0.1;
			TVector3 pos_plane_2 = pos;
			X = pos.X();
			Y = pos.Y();
			RR = X*X+Y*Y;
			if (Rmin_cuts[2] < RR && RR < Rmax_cuts[2]){//(Rmin < RR && RR < Rmax){
				TVector3 pos_reco_plane_2(stepfunction(&X,&pixelsize),
						stepfunction(&Y,&pixelsize),
						11.800);
				int imodule = (int) floor(((pos_reco_plane_2.Phi()+3.1415)/(2*3.1415))*10.);
				if (imodule > 9){ /*cout << imodule << endl;*/ imodule = 0; }
				if (imodule < 0){ cout << imodule << endl; imodule = 0; }
				count_2[imodule]++;
				if (imodule == 0 || imodule == 3 || imodule == 5 || imodule == 8)
					hist_XY_plane_2.Fill(pos_reco_plane_2.X(), pos_reco_plane_2.Y());
			}
			// ****** plane 3 ********
			// smear angle due to silicon in the xz plane
			tan_theta_new = tan(atan(dir.X()/dir.Z())+gRandom->Gaus(0,theta_0_slicon));
			l_x2 = dir.X()*dir.X()+dir.Z()*dir.Z();
			x_new = sqrt(l_x2 / (1/(tan_theta_new*tan_theta_new)+1));
			// recover the lost sign
			if (dir.X() < 0 && x_new > 0) x_new = -x_new;
			z_new = sqrt(l_x2 - x_new*x_new);
			dir.SetX(x_new);
			dir.SetZ(z_new);
			// same smearing procedure in the yz plane
			tan_theta_new = tan(atan(dir.Y()/dir.Z())+gRandom->Gaus(0,theta_0_slicon));
			l_y2 = dir.Y()*dir.Y()+dir.Z()*dir.Z();
			y_new = sqrt(l_y2 / (1/(tan_theta_new*tan_theta_new)+1));
			// recover the lost sign
			if (dir.Y() < 0 && y_new > 0) y_new = -y_new;
			z_new = sqrt(l_y2 - y_new*y_new);
			dir.SetY(y_new);
			dir.SetZ(z_new);
			// propagate to the second plane
			pos += dir.Unit()*0.1;
			TVector3 pos_plane_3 = pos;
			X = pos.X();
			Y = pos.Y();
			RR = X*X+Y*Y;
			if (Rmin_cuts[3] < RR && RR < Rmax_cuts[3]){//(Rmin < RR && RR < Rmax){
				TVector3 pos_reco_plane_3(stepfunction(&X,&pixelsize),
						stepfunction(&Y,&pixelsize),
						11.900);
				int imodule = (int) floor(((pos_reco_plane_3.Phi()+3.1415)/(2*3.1415))*10.);
				if (imodule > 9){ /*cout << imodule << endl;*/ imodule = 0; }
				if (imodule < 0){ cout << imodule << endl; imodule = 0; }
				count_3[imodule]++;
				if (imodule == 0 || imodule == 3 || imodule == 5 || imodule == 8)
					hist_XY_plane_3.Fill(pos_reco_plane_3.X(), pos_reco_plane_3.Y());
			}

		}

		cout << " analyzing histograms " << endl;
		cout << " mean pos of ip was (generated)" << endl;
		cout << " x = " << hist_gen_x_y.GetMean(1)*1e3 << " mm ("<< offsetx <<" mm)" << endl;
		cout << " y = " << hist_gen_x_y.GetMean(2)*1e3 << " mm ("<< offsety <<" mm)" << endl;

		cout << " mean dir of beam was (generated)" << endl;
		cout << " theta = " << hist_gen_theta_phi.GetMean(1)*1e3 << " mrad ("<< tiltx <<" mrad)" << endl;
		cout << " phi = " << hist_gen_theta_phi.GetMean(2) << " rad ("<< tilty <<" mrad)" << endl;
		cout << " dx/dz = " << hist_gen_dx_dy.GetMean(1)*1e3 << " mrad ("<< tiltx <<" mrad)" << endl;
		cout << " dy/dz = " << hist_gen_dx_dy.GetMean(2)*1e3 << " mrad ("<< tilty <<" mrad)" << endl;

		cout << " mean in plane 0,1,2,3 is " << endl;
		vector < double > means_x;
		vector < double > widths_x;
		vector < double > means_y;
		vector < double > widths_y;
		vector < double > z_err; z_err.push_back(0); z_err.push_back(0); z_err.push_back(0); z_err.push_back(0);
		means_x.push_back(hist_XY_plane_0.GetMean(1)); meanx_0 = means_x[0];
		means_x.push_back(hist_XY_plane_1.GetMean(1)); meanx_1 = means_x[1];
		means_x.push_back(hist_XY_plane_2.GetMean(1)); meanx_2 = means_x[2];
		means_x.push_back(hist_XY_plane_3.GetMean(1)); meanx_3 = means_x[3];
		means_y.push_back(hist_XY_plane_0.GetMean(2)); meany_0 = means_x[0];
		means_y.push_back(hist_XY_plane_1.GetMean(2)); meany_1 = means_x[1];
		means_y.push_back(hist_XY_plane_2.GetMean(2)); meany_2 = means_x[2];
		means_y.push_back(hist_XY_plane_3.GetMean(2)); meany_3 = means_x[3];

		widths_x.push_back(hist_XY_plane_0.GetRMS(1)/sqrt(hist_XY_plane_0.GetEntries())); rmsx_0 = hist_XY_plane_0.GetRMS(1);
		widths_x.push_back(hist_XY_plane_1.GetRMS(1)/sqrt(hist_XY_plane_1.GetEntries())); rmsx_1 = hist_XY_plane_1.GetRMS(1);
		widths_x.push_back(hist_XY_plane_2.GetRMS(1)/sqrt(hist_XY_plane_2.GetEntries())); rmsx_2 = hist_XY_plane_2.GetRMS(1);
		widths_x.push_back(hist_XY_plane_3.GetRMS(1)/sqrt(hist_XY_plane_3.GetEntries())); rmsx_3 = hist_XY_plane_3.GetRMS(1);
		widths_y.push_back(hist_XY_plane_0.GetRMS(2)/sqrt(hist_XY_plane_0.GetEntries())); rmsy_0 = hist_XY_plane_0.GetRMS(2);
		widths_y.push_back(hist_XY_plane_1.GetRMS(2)/sqrt(hist_XY_plane_1.GetEntries())); rmsy_1 = hist_XY_plane_1.GetRMS(2);
		widths_y.push_back(hist_XY_plane_2.GetRMS(2)/sqrt(hist_XY_plane_2.GetEntries())); rmsy_2 = hist_XY_plane_2.GetRMS(2);
		widths_y.push_back(hist_XY_plane_3.GetRMS(2)/sqrt(hist_XY_plane_3.GetEntries())); rmsy_3 = hist_XY_plane_3.GetRMS(2);

		//TMarker markeroffset(offsetx, offsety, 4);
		//markeroffset.SetMarkerColor(2);
		canvas_distributions.cd(1);
		hist_XY_plane_0.Draw("COLZ");
		//markeroffset.Draw();
		gPad->Update();
		canvas_distributions.cd(2);
		hist_XY_plane_1.Draw("COLZ");
		//markeroffset.Draw();
		gPad->Update();
		canvas_distributions.cd(3);
		hist_XY_plane_2.Draw("COLZ");
		//markeroffset.Draw();
		gPad->Update();
		canvas_distributions.cd(4);
		hist_XY_plane_3.Draw("COLZ");
		double ip_size = 1.; // mm
		double ip_div = 1.; // mrad
		double ip_width_z = 5.; // mm
		// z positions of planes in m
		vector < double > z; z.push_back(11.5); z.push_back(11.7); z.push_back(11.8); z.push_back(11.9);
		// acceptance in the plane (already the power of 2)
		const double Rmin = 0.035*0.035;
		const double Rmax = (0.035+0.06)*(0.035+0.06);
		// cut for a cone starting from ip
		// angles for a cone cut
		const double theta_cut_low = 4.; // mrad
		const double theta_cut_high = 7.; // mrad
		vector < double > Rmin_cuts;
		vector < double > Rmax_cuts;
		for (int iplane = 0; iplane < 4; iplane++ ){
			Rmin_cuts.push_back(theta_cut_low*1e-3*z[iplane]*theta_cut_low*1e-3*z[iplane]);
			Rmax_cuts.push_back(theta_cut_high*1e-3*z[iplane]*theta_cut_high*1e-3*z[iplane]);
		}

		TCanvas canvas("canvas", "canvas", 600, 1000);
		canvas.Divide(2,1);
		canvas.cd(1);
		TCanvas canvas_distributions("canvas_distributions", "distributions", 800, 800);
		canvas_distributions.Divide(2,2);

		//TFile treefile("asymmetries_results.root", "RECREATE");
		TTree results("asymmetry_results", "Asymmetry results");
		double offsetx, offsety, tiltx, tilty, momentum,
			offsetx_fit, offsety_fit, tiltx_fit, tilty_fit,
			meanx_0, meany_0, rmsx_0, rmsy_0,
			meanx_1, meany_1, rmsx_1, rmsy_1,
			meanx_2, meany_2, rmsx_2, rmsy_2,
			meanx_3, meany_3, rmsx_3, rmsy_3;
		int	count_0[10] = {0,0,0,0,0,0,0,0,0,0};
		int count_1[10] = {0,0,0,0,0,0,0,0,0,0};
		int count_2[10] = {0,0,0,0,0,0,0,0,0,0};
		int count_3[10] = {0,0,0,0,0,0,0,0,0,0};
		TLine line;

		results.Branch("offsetx", &offsetx);
		results.Branch("offsety", &offsety);
		results.Branch("tiltx", &tiltx);
		results.Branch("tilty", &tilty);
		results.Branch("momentum", &momentum);
		results.Branch("offsetx_fit", &offsetx_fit);
		results.Branch("offsety_fit", &offsety_fit);
		results.Branch("tiltx_fit", &tiltx_fit);
		results.Branch("tilty_fit", &tilty_fit);
		results.Branch("meanx_0", &meanx_0);
		results.Branch("meany_0", &meany_0);
		results.Branch("rmsx_0", &rmsx_0);
		results.Branch("rmsy_0", &rmsy_0);
		results.Branch("meanx_1", &meanx_1);
		results.Branch("meany_1", &meany_1);
		results.Branch("rmsx_1", &rmsx_1);
		results.Branch("rmsy_1", &rmsy_1);
		results.Branch("meanx_2", &meanx_2);
		results.Branch("meany_2", &meany_2);
		results.Branch("rmsx_2", &rmsx_2);
		results.Branch("rmsy_2", &rmsy_2);
		results.Branch("meanx_3", &meanx_3);
		results.Branch("meany_3", &meany_3);
		results.Branch("rmsx_3", &rmsx_3);
		results.Branch("rmsy_3", &rmsy_3);
		results.Branch("count_0", count_0, "count_0[10]/I");
		results.Branch("count_1", count_1, "count_1[10]/I");
		results.Branch("count_2", count_2, "count_2[10]/I");
		results.Branch("count_3", count_3, "count_3[10]/I");

		// calculate the number of loops
		int nloops = 0;
		for (offsetx = -2.; offsetx <= 2.; offsetx+=2.) // mm
		for (offsety = -2.; offsety <= 2.; offsety+=2.)
		for (tiltx = -.5; tiltx <= .5; tiltx+=.5) // mrad
		for (tilty = -.5; tilty <= .5; tilty+=.5) // mrad
		for (int imom = 0; imom < nmomenta; imom++){
			nloops++;
		}

		int iloop = 0;

		TGraphAsymmErrors* graph_model = NULL;
		TH1* model_distribution = NULL;
		for (int imom = 0; imom < nmomenta; imom++){
			canvas.cd(1);
			momentum = momenta[imom];
			// try to retrieve the DPM density function to model the angular distribution
			stringstream modelfilename;
			modelfilename << "DPMModels_" << momentum << "_model.root";
			TFile model(modelfilename.str().c_str(), "READ");
			delete graph_model;
			TGraphAsymmErrors* graph_model = (TGraphAsymmErrors*) model.Get("Graph");
			TH1F model_distribution("model_distribution", "distribution", 1000, 1., 10.);
			model.Close();
			if (!graph_model){
				cout << " Warning: no density function found in " << modelfilename.str() << endl;
			} else {
				//graph_model->Draw("ALP");
				//gPad->Update();
				for (double ibin = 1.; ibin <= model_distribution.GetXaxis()->GetNbins(); ibin++){
					double theta = model_distribution.GetXaxis()->GetBinCenter(ibin); // mrad
					model_distribution.Fill(theta, graph_model->Eval(theta));
				}
				gPad->Clear();
				model_distribution.Draw();
				gPad->Update();
			}
		for (offsetx = -2.; offsetx <= 2.; offsetx+=2.){ // mm
		for (offsety = -2.; offsety <= 2.; offsety+=2.){
		for (tiltx = -.5; tiltx <= .5; tiltx+=.5){ // mrad
		for (tilty = -.5; tilty <= .5; tilty+=.5){ // mrad
			iloop++;
			cout << " **** " << iloop << " of " << nloops << " ******** " << endl;
			cout << " **** " << iloop*100/nloops << "% **********" << endl;
			stringstream filename;
			filename << "asymmetry_studies_"<< momentum << "_GeVc_x" << offsetx << "_mm_y_" << offsety << "_mm_dxdz_" << tiltx << "_mrad_dydz_" << tilty << "_mrad";
			//TFile file((filename.str()+".root").c_str(), "RECREATE");

			TH2F hist_gen_theta_phi ("hist_gen_theta_phi", "theta phi generated distribution", 100, 2e-3, 9e-3, 100, -3.141, 3.141);
			TH2F hist_gen_x_y ("hist_gen_x_y", "x y generated distribution", 100, -5e-3, 5e-3, 100, -5e-3, 5e-3);
			TH2F hist_gen_dx_dy ("hist_gen_dx_dy", "dx/dz dy/dz generated distribution", 100, -15e-3, 15e-3, 100, -15e-3, 15e-3);
			TH2F hist_XY_plane_0 ("hist_XY_plane_0", "xy distribution at plane 0", 100, -0.100, 0.100, 100, -0.100, 0.100);
			TH2F hist_XY_plane_1 ("hist_XY_plane_1", "xy distribution at plane 1", 100, -0.100, 0.100, 100, -0.100, 0.100);
			TH2F hist_XY_plane_2 ("hist_XY_plane_2", "xy distribution at plane 2", 100, -0.100, 0.100, 100, -0.100, 0.100);
			TH2F hist_XY_plane_3 ("hist_XY_plane_3", "xy distribution at plane 3", 100, -0.100, 0.100, 100, -0.100, 0.100);
			//TH2F hist_phi_count("hist_phi_count", "counts over plane number", 4, -0.5, 3.5, 10, -3.141, 3.141);
			//TH2F hist_center_of_gravity("hist_center_of_gravity", "center of gravity over plane number", 4, -0.5, 3.5, 10, -3.141, 3.141);

			// calculate some multiple scattering constants

			double beta = 1/sqrt(1+pow(0.938/momentum,2));
			cout << "the beta factor for " << momentum << " GeV/c antiprotons is " << beta << endl;
			double x_over_x0 = 2.329/21.82*350/1000/10;
			cout << "x/X0 is" << x_over_x0 << endl;
			double theta_0_slicon = 13.6/beta/(momentum*1e3)*sqrt(x_over_x0)*(1+0.038*log(x_over_x0));
			cout << "scattering angle of 350 mum silicon is " << theta_0_slicon << endl;
			// assuming a multiple scattering effect in the transition region of a silicon plane of half thickness
			double theta_0_kapton = 13.6/beta/(momentum*1e3)*sqrt(x_over_x0/2.)*(1+0.038*log(x_over_x0/2.));
			for (int ievent = 0; ievent < 1000000*1.7; ievent++){
				double phi =  gRandom->Uniform(-3.141, 3.141);
				double theta;
				if (graph_model) theta = model_distribution.GetRandom()*1e-3;
				else theta = gRandom->Uniform(1e-3, 10e-3);
				TVector3 dir(0., 0., momentum);
				dir.SetTheta(theta);
				dir.SetPhi(phi);
				TVector3 pos(gRandom->Gaus(offsetx, ip_size)*1e-3,gRandom->Gaus(offsety, ip_size)*1e-3,gRandom->Uniform(-ip_width_z, ip_width_z)*1e-3); // m
				// rotate according to a tilted and smeared beam
				TVector3 dir_beam(gRandom->Gaus(tiltx, ip_div)*1e-3,gRandom->Gaus(tilty, ip_div)*1e-3,1.);
				//cout << "dir_beam " << dir_beam.X()/dir_beam.Z()-tiltx*1e-3 << " " << dir_beam.Y()/dir_beam.Z()-tilty*1e-3 << endl;
				dir.RotateUz(dir_beam.Unit());
				//cout << "dir " << (dir.X()/dir.Z()*dir.X()/dir.Z()+dir.Y()/dir.Z()*dir.Y()/dir.Z()) tiltx*1e-3 << " " << dir.Y()/dir.Z()-tilty*1e-3 << endl;
				hist_gen_x_y.Fill(pos.X(), pos.Y());
				hist_gen_theta_phi.Fill(dir.Theta(), dir.Phi());
				hist_gen_dx_dy.Fill(dir.X()/dir.Z(), dir.Y()/dir.Z());

				//hist_gen.Fill(dir.Theta(), dir.Phi());
				// propagate to the transition region at about 11m
				pos += dir.Unit()*11.; // sure propagation is not fully correct as it is not propagated to a plane
				// smear angle due to silicon in the xz plane
				double tan_theta_new = tan(atan(dir.X()/dir.Z())+gRandom->Gaus(0,theta_0_kapton));
				double l_x2 = dir.X()*dir.X()+dir.Z()*dir.Z();
				double x_new = sqrt(l_x2 / (1/(tan_theta_new*tan_theta_new)+1));
				double z_new = sqrt(l_x2 - x_new*x_new);
				// recover the lost sign
				if (dir.X() < 0 && x_new > 0) x_new = -x_new;
				dir.SetX(x_new);
				dir.SetZ(z_new);
				// same smearing procedure in the yz plane
				tan_theta_new = tan(atan(dir.Y()/dir.Z())+gRandom->Gaus(0,theta_0_kapton));
				double l_y2 = dir.Y()*dir.Y()+dir.Z()*dir.Z();
				double y_new = sqrt(l_y2 / (1/(tan_theta_new*tan_theta_new)+1));
				// recover the lost sign
				if (dir.Y() < 0 && y_new > 0) y_new = -y_new;
				z_new = sqrt(l_y2 - y_new*y_new);
				dir.SetY(y_new);
				dir.SetZ(z_new);
				// propagate to the first plane
				pos += dir.Unit()*0.5;
				TVector3 pos_plane_0 = pos;
				//hist_phi_count.Fill(0, pos_plane_0.Phi());
				// simulate a granularity of the detector in theta of XX mum
				double X = pos.X();
				double Y = pos.Y();
				double RR = X*X+Y*Y;
				if (Rmin_cuts[0] < RR && RR < Rmax_cuts[0]){//(Rmin < RR && RR < Rmax){
					TVector3 pos_reco_plane_0(stepfunction(&X,&pixelsize),
							stepfunction(&Y,&pixelsize),
							11.500);
					int imodule = (int) floor(((pos_reco_plane_0.Phi()+3.1415)/(2*3.1415))*10.);
					if (imodule > 9){ /*cout << imodule << endl;*/ imodule = 0; }
					if (imodule < 0){ cout << imodule << endl; imodule = 0; }
					count_0[imodule]++;
					if (imodule == 0 || imodule == 3 || imodule == 5 || imodule == 8)
					hist_XY_plane_0.Fill(pos_reco_plane_0.X(), pos_reco_plane_0.Y());
				}
				// ****** plane 1 ********
				// smear angle due to silicon in the xz plane
				tan_theta_new = tan(atan(dir.X()/dir.Z())+gRandom->Gaus(0,theta_0_slicon));
				l_x2 = dir.X()*dir.X()+dir.Z()*dir.Z();
				x_new = sqrt(l_x2 / (1/(tan_theta_new*tan_theta_new)+1));
				// recover the lost sign
				if (dir.X() < 0 && x_new > 0) x_new = -x_new;
				z_new = sqrt(l_x2 - x_new*x_new);
				dir.SetX(x_new);
				dir.SetZ(z_new);
				// same smearing procedure in the yz plane
				tan_theta_new = tan(atan(dir.Y()/dir.Z())+gRandom->Gaus(0,theta_0_slicon));
				l_y2 = dir.Y()*dir.Y()+dir.Z()*dir.Z();
				y_new = sqrt(l_y2 / (1/(tan_theta_new*tan_theta_new)+1));
				// recover the lost sign
				if (dir.Y() < 0 && y_new > 0) y_new = -y_new;
				z_new = sqrt(l_y2 - y_new*y_new);
				dir.SetY(y_new);
				dir.SetZ(z_new);
				// propagate to the second plane
				pos += dir.Unit()*0.2;
				TVector3 pos_plane_1 = pos;
				X = pos.X();
				Y = pos.Y();
				RR = X*X+Y*Y;
				if (Rmin_cuts[1] < RR && RR < Rmax_cuts[1]){//(Rmin < RR && RR < Rmax){
					TVector3 pos_reco_plane_1(stepfunction(&X,&pixelsize),
							stepfunction(&Y,&pixelsize),
							11.700);
					int imodule = (int) floor(((pos_reco_plane_1.Phi()+3.1415)/(2*3.1415))*10.);
					if (imodule > 9){ /*cout << imodule << endl;*/ imodule = 0; }
					if (imodule < 0){ cout << imodule << endl; imodule = 0; }
					count_1[imodule]++;
					if (imodule == 0 || imodule == 3 || imodule == 5 || imodule == 8)
						hist_XY_plane_1.Fill(pos_reco_plane_1.X(), pos_reco_plane_1.Y());
				}
				// ****** plane 2 ********
				// smear angle due to silicon in the xz plane
				tan_theta_new = tan(atan(dir.X()/dir.Z())+gRandom->Gaus(0,theta_0_slicon));
				l_x2 = dir.X()*dir.X()+dir.Z()*dir.Z();
				x_new = sqrt(l_x2 / (1/(tan_theta_new*tan_theta_new)+1));
				// recover the lost sign
				if (dir.X() < 0 && x_new > 0) x_new = -x_new;
				z_new = sqrt(l_x2 - x_new*x_new);
				dir.SetX(x_new);
				dir.SetZ(z_new);
				// same smearing procedure in the yz plane
				tan_theta_new = tan(atan(dir.Y()/dir.Z())+gRandom->Gaus(0,theta_0_slicon));
				l_y2 = dir.Y()*dir.Y()+dir.Z()*dir.Z();
				y_new = sqrt(l_y2 / (1/(tan_theta_new*tan_theta_new)+1));
				// recover the lost sign
				if (dir.Y() < 0 && y_new > 0) y_new = -y_new;
				z_new = sqrt(l_y2 - y_new*y_new);
				dir.SetY(y_new);
				dir.SetZ(z_new);
				// propagate to the second plane
				pos += dir.Unit()*0.1;
				TVector3 pos_plane_2 = pos;
				X = pos.X();
				Y = pos.Y();
				RR = X*X+Y*Y;
				if (Rmin_cuts[2] < RR && RR < Rmax_cuts[2]){//(Rmin < RR && RR < Rmax){
					TVector3 pos_reco_plane_2(stepfunction(&X,&pixelsize),
							stepfunction(&Y,&pixelsize),
							11.800);
					int imodule = (int) floor(((pos_reco_plane_2.Phi()+3.1415)/(2*3.1415))*10.);
					if (imodule > 9){ /*cout << imodule << endl;*/ imodule = 0; }
					if (imodule < 0){ cout << imodule << endl; imodule = 0; }
					count_2[imodule]++;
					if (imodule == 0 || imodule == 3 || imodule == 5 || imodule == 8)
						hist_XY_plane_2.Fill(pos_reco_plane_2.X(), pos_reco_plane_2.Y());
				}
				// ****** plane 3 ********
				// smear angle due to silicon in the xz plane
				tan_theta_new = tan(atan(dir.X()/dir.Z())+gRandom->Gaus(0,theta_0_slicon));
				l_x2 = dir.X()*dir.X()+dir.Z()*dir.Z();
				x_new = sqrt(l_x2 / (1/(tan_theta_new*tan_theta_new)+1));
				// recover the lost sign
				if (dir.X() < 0 && x_new > 0) x_new = -x_new;
				z_new = sqrt(l_x2 - x_new*x_new);
				dir.SetX(x_new);
				dir.SetZ(z_new);
				// same smearing procedure in the yz plane
				tan_theta_new = tan(atan(dir.Y()/dir.Z())+gRandom->Gaus(0,theta_0_slicon));
				l_y2 = dir.Y()*dir.Y()+dir.Z()*dir.Z();
				y_new = sqrt(l_y2 / (1/(tan_theta_new*tan_theta_new)+1));
				// recover the lost sign
				if (dir.Y() < 0 && y_new > 0) y_new = -y_new;
				z_new = sqrt(l_y2 - y_new*y_new);
				dir.SetY(y_new);
				dir.SetZ(z_new);
				// propagate to the second plane
				pos += dir.Unit()*0.1;
				TVector3 pos_plane_3 = pos;
				X = pos.X();
				Y = pos.Y();
				RR = X*X+Y*Y;
				if (Rmin_cuts[3] < RR && RR < Rmax_cuts[3]){//(Rmin < RR && RR < Rmax){
					TVector3 pos_reco_plane_3(stepfunction(&X,&pixelsize),
							stepfunction(&Y,&pixelsize),
							11.900);
					int imodule = (int) floor(((pos_reco_plane_3.Phi()+3.1415)/(2*3.1415))*10.);
					if (imodule > 9){ /*cout << imodule << endl;*/ imodule = 0; }
					if (imodule < 0){ cout << imodule << endl; imodule = 0; }
					count_3[imodule]++;
					if (imodule == 0 || imodule == 3 || imodule == 5 || imodule == 8)
						hist_XY_plane_3.Fill(pos_reco_plane_3.X(), pos_reco_plane_3.Y());
				}

			}

			cout << " analyzing histograms " << endl;
			cout << " mean pos of ip was (generated)" << endl;
			cout << " x = " << hist_gen_x_y.GetMean(1)*1e3 << " mm ("<< offsetx <<" mm)" << endl;
			cout << " y = " << hist_gen_x_y.GetMean(2)*1e3 << " mm ("<< offsety <<" mm)" << endl;

			cout << " mean dir of beam was (generated)" << endl;
			cout << " theta = " << hist_gen_theta_phi.GetMean(1)*1e3 << " mrad ("<< tiltx <<" mrad)" << endl;
			cout << " phi = " << hist_gen_theta_phi.GetMean(2) << " rad ("<< tilty <<" mrad)" << endl;
			cout << " dx/dz = " << hist_gen_dx_dy.GetMean(1)*1e3 << " mrad ("<< tiltx <<" mrad)" << endl;
			cout << " dy/dz = " << hist_gen_dx_dy.GetMean(2)*1e3 << " mrad ("<< tilty <<" mrad)" << endl;

			cout << " mean in plane 0,1,2,3 is " << endl;
			vector < double > means_x;
			vector < double > widths_x;
			vector < double > means_y;
			vector < double > widths_y;
			vector < double > z_err; z_err.push_back(0); z_err.push_back(0); z_err.push_back(0); z_err.push_back(0);
			means_x.push_back(hist_XY_plane_0.GetMean(1)); meanx_0 = means_x[0];
			means_x.push_back(hist_XY_plane_1.GetMean(1)); meanx_1 = means_x[1];
			means_x.push_back(hist_XY_plane_2.GetMean(1)); meanx_2 = means_x[2];
			means_x.push_back(hist_XY_plane_3.GetMean(1)); meanx_3 = means_x[3];
			means_y.push_back(hist_XY_plane_0.GetMean(2)); meany_0 = means_x[0];
			means_y.push_back(hist_XY_plane_1.GetMean(2)); meany_1 = means_x[1];
			means_y.push_back(hist_XY_plane_2.GetMean(2)); meany_2 = means_x[2];
			means_y.push_back(hist_XY_plane_3.GetMean(2)); meany_3 = means_x[3];

			widths_x.push_back(hist_XY_plane_0.GetRMS(1)/sqrt(hist_XY_plane_0.GetEntries())); rmsx_0 = hist_XY_plane_0.GetRMS(1);
			widths_x.push_back(hist_XY_plane_1.GetRMS(1)/sqrt(hist_XY_plane_1.GetEntries())); rmsx_1 = hist_XY_plane_1.GetRMS(1);
			widths_x.push_back(hist_XY_plane_2.GetRMS(1)/sqrt(hist_XY_plane_2.GetEntries())); rmsx_2 = hist_XY_plane_2.GetRMS(1);
			widths_x.push_back(hist_XY_plane_3.GetRMS(1)/sqrt(hist_XY_plane_3.GetEntries())); rmsx_3 = hist_XY_plane_3.GetRMS(1);
			widths_y.push_back(hist_XY_plane_0.GetRMS(2)/sqrt(hist_XY_plane_0.GetEntries())); rmsy_0 = hist_XY_plane_0.GetRMS(2);
			widths_y.push_back(hist_XY_plane_1.GetRMS(2)/sqrt(hist_XY_plane_1.GetEntries())); rmsy_1 = hist_XY_plane_1.GetRMS(2);
			widths_y.push_back(hist_XY_plane_2.GetRMS(2)/sqrt(hist_XY_plane_2.GetEntries())); rmsy_2 = hist_XY_plane_2.GetRMS(2);
			widths_y.push_back(hist_XY_plane_3.GetRMS(2)/sqrt(hist_XY_plane_3.GetEntries())); rmsy_3 = hist_XY_plane_3.GetRMS(2);

			//TMarker markeroffset(offsetx, offsety, 4);
			//markeroffset.SetMarkerColor(2);
			canvas_distributions.cd(1);
			hist_XY_plane_0.Draw("COLZ");
			//markeroffset.Draw();
			gPad->Update();
			canvas_distributions.cd(2);
			hist_XY_plane_1.Draw("COLZ");
			//markeroffset.Draw();
			gPad->Update();
			canvas_distributions.cd(3);
			hist_XY_plane_2.Draw("COLZ");
			//markeroffset.Draw();
			gPad->Update();
			canvas_distributions.cd(4);
			hist_XY_plane_3.Draw("COLZ");
			//markeroffset.Draw();
			gPad->Update();

			// save values to tree


			cout << " x = " << meanx_0 << ", "<< meanx_1 << ", "<< meanx_2 << ", "<< meanx_3 << endl;
			cout << " y = " << meany_0 << ", "<< meany_1 << ", "<< meany_2 << ", "<< meany_3 << endl;
			cout << " rms x = " << rmsx_0 << ", "<< rmsx_1 << ", "<< rmsx_2 << ", "<< rmsx_3 << endl;
			cout << " rms y = " << rmsy_0 << ", "<< rmsy_1 << ", "<< rmsy_2 << ", "<< rmsy_3 << endl;

			TGraphErrors graph_average_x(means_x.size(), &z[0], &means_x[0], NULL, &widths_x[0]);
			graph_average_x.SetNameTitle("graph_average_x", "average x over z");
			canvas.cd(1);
			gPad->Clear();
			graph_average_x.Draw("ALP");
			graph_average_x.GetYaxis()->SetRangeUser(-0.02, 0.02);
			TF1 funcx("funcx", "pol1", -0.1, 0.5);
			graph_average_x.Fit(&funcx);
			line.DrawLine(z[0], offsetx*1e-3, z[3], offsetx*1e-3);
			line.DrawLine(z[0], tiltx*1e-3*z[0]+offsetx*1e-3, z[3], tiltx*1e-3*z[3]+offsetx*1e-3);
			gPad->Update();
			//sleep(2);
			//graph_average_x.Write();
			canvas.cd(2);
			TGraphErrors graph_average_y(means_y.size(), &z[0], &means_y[0], NULL, &widths_y[0]);
			graph_average_y.SetNameTitle("graph_average_y", "average y over z");
			gPad->Clear();
			graph_average_y.Draw("ALP");
			graph_average_y.GetYaxis()->SetRangeUser(-0.02, 0.02);
			TF1 funcy("funcy", "pol1", -0.1, 0.5);
			graph_average_y.Fit(&funcy);
			line.DrawLine(z[0], offsety*1e-3, z[3], offsety*1e-3);
			line.DrawLine(z[0], tilty*1e-3*z[0]+offsety*1e-3, z[3], tilty*1e-3*z[3]+offsety*1e-3);
			//graph_average_y.Write();
			gPad->Update();

			cout << " mean pos of ip fit (generated)" << endl;
			offsetx_fit = funcx.GetParameter(0)*1e3;
			offsety_fit = funcy.GetParameter(0)*1e3;
			cout << " x = " << offsetx_fit << " mm ("<< offsetx <<" mm)" << endl;
			cout << " y = " << offsety_fit << " mm ("<< offsety <<" mm)" << endl;

			tiltx_fit = funcx.GetParameter(1)*1e3;
			tilty_fit = funcy.GetParameter(1)*1e3;
			cout << " mean dir of beam fit (generated)" << endl;
			cout << " dx/dz = " << tiltx_fit << " mrad ("<< tiltx <<" mrad)" << endl;
			cout << " dy/dz = " << tilty_fit << " mrad ("<< tilty <<" mrad)" << endl;

			cout << endl;
			results.Fill();
			//file.Write();
			//file.Close();
		}
		}
		}
		}
		}
		TFile treefile("asymmetries_results.root", "RECREATE");
		results.SetDirectory(&treefile);
		results.Write();
		treefile.Close();
		cout << " done " << endl;
		//markeroffset.Draw();
		gPad->Update();

		// save values to tree


		cout << " x = " << meanx_0 << ", "<< meanx_1 << ", "<< meanx_2 << ", "<< meanx_3 << endl;
		cout << " y = " << meany_0 << ", "<< meany_1 << ", "<< meany_2 << ", "<< meany_3 << endl;
		cout << " rms x = " << rmsx_0 << ", "<< rmsx_1 << ", "<< rmsx_2 << ", "<< rmsx_3 << endl;
		cout << " rms y = " << rmsy_0 << ", "<< rmsy_1 << ", "<< rmsy_2 << ", "<< rmsy_3 << endl;

		TGraphErrors graph_average_x(means_x.size(), &z[0], &means_x[0], NULL, &widths_x[0]);
		graph_average_x.SetNameTitle("graph_average_x", "average x over z");
		canvas.cd(1);
		gPad->Clear();
		graph_average_x.Draw("ALP");
		graph_average_x.GetYaxis()->SetRangeUser(-0.02, 0.02);
		TF1 funcx("funcx", "pol1", -0.1, 0.5);
		graph_average_x.Fit(&funcx);
		line.DrawLine(z[0], offsetx*1e-3, z[3], offsetx*1e-3);
		line.DrawLine(z[0], tiltx*1e-3*z[0]+offsetx*1e-3, z[3], tiltx*1e-3*z[3]+offsetx*1e-3);
		gPad->Update();
		//sleep(2);
		//graph_average_x.Write();
		canvas.cd(2);
		TGraphErrors graph_average_y(means_y.size(), &z[0], &means_y[0], NULL, &widths_y[0]);
		graph_average_y.SetNameTitle("graph_average_y", "average y over z");
		gPad->Clear();
		graph_average_y.Draw("ALP");
		graph_average_y.GetYaxis()->SetRangeUser(-0.02, 0.02);
		TF1 funcy("funcy", "pol1", -0.1, 0.5);
		graph_average_y.Fit(&funcy);
		line.DrawLine(z[0], offsety*1e-3, z[3], offsety*1e-3);
		line.DrawLine(z[0], tilty*1e-3*z[0]+offsety*1e-3, z[3], tilty*1e-3*z[3]+offsety*1e-3);
		//graph_average_y.Write();
		gPad->Update();

		cout << " mean pos of ip fit (generated)" << endl;
		offsetx_fit = funcx.GetParameter(0)*1e3;
		offsety_fit = funcy.GetParameter(0)*1e3;
		cout << " x = " << offsetx_fit << " mm ("<< offsetx <<" mm)" << endl;
		cout << " y = " << offsety_fit << " mm ("<< offsety <<" mm)" << endl;

		tiltx_fit = funcx.GetParameter(1)*1e3;
		tilty_fit = funcy.GetParameter(1)*1e3;
		cout << " mean dir of beam fit (generated)" << endl;
		cout << " dx/dz = " << tiltx_fit << " mrad ("<< tiltx <<" mrad)" << endl;
		cout << " dy/dz = " << tilty_fit << " mrad ("<< tilty <<" mrad)" << endl;

		cout << endl;
		results.Fill();
		//file.Write();
		//file.Close();
	}
	}
	}
	}
	}
	TFile treefile("asymmetries_results.root", "RECREATE");
	results.SetDirectory(&treefile);
	results.Write();
	treefile.Close();
	cout << " done " << endl;
}


void MickeyMouse_MC_resolution(){
	const int nmomenta = 5;
	double momentum[nmomenta] = {1.5, 4.06, 8.9, 11.91, 15.00};
	double pixelsize = 80;

	for (int imomentum = 0; imomentum < nmomenta; imomentum++){

		stringstream filename;
		filename << "resolution_studies_"<<pixelsize<<"mu_" << momentum[imomentum] << "_GeV";
		TFile file((filename.str()+".root").c_str(), "RECREATE");
		TH2F hist_gen ("hist_gen", "generated distribution", 100, 2e-3, 9e-3, 100, -3.141, 3.141);
		TH2F hist_XY_plane_0 ("hist_XY_plane_0", "xy distribution at plane 0", 100, -0.100, 0.100, 100, -0.100, 0.100);
		TH2F hist_XY_plane_1 ("hist_XY_plane_1", "xy distribution at plane 1", 100, -0.100, 0.100, 100, -0.100, 0.100);
		TH1F hist_theta_reso("hist_theta_reso", "intrinsic #Theta resolution", 1000, -1e-3, 1e-3);
		TH1F hist_theta_reso_reco("hist_theta_reso_reco", "pixelized #Theta resolution", 1000, -1e-3, 1e-3);
		TH2F hist_ana_crosscheck("hist_ana_crosscheck", "anastasia crosscheck", 100, -400, 400, 100, -15, 15);
		//TF1* testfunction = new TF1("testfunction", &stepfunction, -1.e-3, 1.e-3, 1);
		//testfunction->SetNpx(10000);
		//testfunction->SetParameter(0,pixelsize);
		//testfunction->Draw("");
		//return;

		TCanvas canvas("canvas", "canvas", 800, 400);
		canvas.Divide(2,1);
		canvas.cd(1);

		// calculate some multiple scattering constants

		double beta = 1/sqrt(1+pow(0.938/momentum[imomentum],2));
		cout << "the beta factor is" << beta << endl;
		double x_over_x0 = 2.329/21.82*150/1000/10;
		cout << "x/X0 is" << x_over_x0 << endl;
		double theta_0_slicon = 13.6/beta/(momentum[imomentum]*1e3)*sqrt(x_over_x0)*(1+0.038*log(x_over_x0));
		cout << "scattering angle of 150 mum silicon is " << theta_0_slicon << endl;
		// assuming a multiple scattering effect in the transition region of a silicon plane of half thickness
		double theta_0_kapton = 13.6/beta/(momentum[imomentum]*1e3)*sqrt(x_over_x0/2.)*(1+0.038*log(x_over_x0/2.));

		for (int ievent = 0; ievent < 1000000; ievent++){
			double phi = gRandom->Uniform(-3.141, 3.141);
			double theta = gRandom->Uniform(3e-3, 8e-3);
			TVector3 pos(0.,0.,0.);
			TVector3 dir(0.,0.,momentum[imomentum]);
			dir.SetTheta(theta);
			dir.SetPhi(phi);
			hist_gen.Fill(dir.Theta(), dir.Phi());
			// propagate to the transition region at about 10m
			pos += dir.Unit()*10; // sure propagation is not fully correct as it is not propagated to a plane
			// smear angle due to silicon in the xz plane
			double tan_theta_new = tan(atan(dir.X()/dir.Z())+gRandom->Gaus(0,theta_0_kapton));
			double l_x2 = dir.X()*dir.X()+dir.Z()*dir.Z();
			double x_new = sqrt(l_x2 / (1/(tan_theta_new*tan_theta_new)+1));
			double z_new = sqrt(l_x2 - x_new*x_new);
			// recover the lost sign
			if (dir.X() < 0 && x_new > 0) x_new = -x_new;
			dir.SetX(x_new);
			dir.SetZ(z_new);
			// same smearing procedure in the yz plane
			tan_theta_new = tan(atan(dir.Y()/dir.Z())+gRandom->Gaus(0,theta_0_kapton));
			double l_y2 = dir.Y()*dir.Y()+dir.Z()*dir.Z();
			double y_new = sqrt(l_y2 / (1/(tan_theta_new*tan_theta_new)+1));
			// recover the lost sign
			if (dir.Y() < 0 && y_new > 0) y_new = -y_new;
			z_new = sqrt(l_y2 - y_new*y_new);
			dir.SetY(y_new);
			dir.SetZ(z_new);
			// propagate to the first plane
			pos += dir.Unit()*0.5;
			TVector3 pos_plane_0 = pos;
			hist_XY_plane_0.Fill(pos_plane_0.X(), pos_plane_0.Y());
			// simulate a granularity of the detector in theta of XX mum
			double X = pos.X();
			double Y = pos.Y();
			TVector3 pos_reco_plane_0(stepfunction(&X,&pixelsize),
					stepfunction(&Y,&pixelsize),
					10.500);
			// smear angle due to silicon in the xz plane
			tan_theta_new = tan(atan(dir.X()/dir.Z())+gRandom->Gaus(0,theta_0_slicon));
			l_x2 = dir.X()*dir.X()+dir.Z()*dir.Z();
			x_new = sqrt(l_x2 / (1/(tan_theta_new*tan_theta_new)+1));
			// recover the lost sign
			if (dir.X() < 0 && x_new > 0) x_new = -x_new;
			z_new = sqrt(l_x2 - x_new*x_new);
			dir.SetX(x_new);
			dir.SetZ(z_new);
			// same smearing procedure in the yz plane
			tan_theta_new = tan(atan(dir.Y()/dir.Z())+gRandom->Gaus(0,theta_0_slicon));
			l_y2 = dir.Y()*dir.Y()+dir.Z()*dir.Z();
			y_new = sqrt(l_y2 / (1/(tan_theta_new*tan_theta_new)+1));
			// recover the lost sign
			if (dir.Y() < 0 && y_new > 0) y_new = -y_new;
			z_new = sqrt(l_y2 - y_new*y_new);
			dir.SetY(y_new);
			dir.SetZ(z_new);
			// propagate to the second plane
			pos += dir.Unit()*0.2;
			TVector3 pos_plane_1 = pos;
			hist_XY_plane_1.Fill(pos_plane_1.X(), pos_plane_1.Y());
			X = pos.X();
			Y = pos.Y();
			TVector3 pos_reco_plane_1(stepfunction(&X,&pixelsize),
					stepfunction(&Y,&pixelsize),
					10.700);
			double theta_propa = (pos_plane_1-pos_plane_0).Theta();
			double theta_reco  = (pos_reco_plane_1-pos_reco_plane_0).Theta();
			double reso_propa = theta_propa-theta;
			hist_theta_reso.Fill(reso_propa);
			double reso_reco = theta_reco-theta;
			hist_theta_reso_reco.Fill(reso_reco);

			TVector3 dir_reco = pos_reco_plane_1-pos_reco_plane_0;
			dir_reco.RotateY(40.0e-3);///3.1416*180.);
			double theta_tilde = sqrt(dir_reco.X()*dir_reco.X()+dir_reco.Y()*dir_reco.Y())/dir_reco.Z()-0.040;
			double phi_tilde   = dir_reco.Y()/dir_reco.X();
			//cout << theta_tilde << " " << phi_tilde << endl;
			hist_ana_crosscheck.Fill(phi_tilde*1e3, theta_tilde*1e3);
		}
		canvas.cd(1);
		hist_theta_reso.Draw();
		gPad->Update();
		canvas.cd(2);
		hist_theta_reso_reco.Draw();
		gPad->Update();
		//hist_ana_crosscheck.Draw();
		//gPad->Update();
		canvas.Print((filename.str()+".pdf").c_str());
		file.Write();
		file.Close();
	}
}

// #include "testclass.h"

// testclass get_testclass(){
// 	//&testclass result = new testclass();
// 	return testclass();


// }

#include "TMatrixD.h"

TMatrixD Get_rotation_matrix(){
	TMatrixD result(3,3);
	result[0][0] = 2.;
	result[1][0] = -1.5;
	result[2][0] = 3.;
	result[0][1] = 6.;
	result[1][1] = -2.5;
	result[2][1] = 2.5;
	result[0][2] = 10.;
	result[1][2] = -0.5;
	result[2][2] = 1.;
	return result;
}

TMatrixD Get_matrix_for_rotation(){
	TMatrixD result(3,3);
	result[0][0] = 1.;
	result[1][0] = 2.5;
	result[2][0] = -2.;
	result[0][1] = -1.;
	result[1][1] = -3.5;
	result[2][1] = 0.5;
	result[0][2] = 4.;
	result[1][2] = 0.5;
	result[2][2] = -4.;
	return result;
}

TMatrixD& Transform_global_to_lmd_local_ref(const TMatrixD& matrix){
	TMatrixD rotmatrix(Get_rotation_matrix());
	return *(new TMatrixD(rotmatrix*matrix*rotmatrix.T()));
}

TMatrixD Transform_global_to_lmd_local_val(const TMatrixD& matrix){
	TMatrixD rotmatrix(Get_rotation_matrix());
	return TMatrixD(rotmatrix*matrix*rotmatrix.T());
}

TMatrixD Transform_global_to_lmd_local_ana(const TMatrixD& matrix){
	TMatrixD rotmatrix(Get_rotation_matrix());
	return TMatrixD(rotmatrix*TMatrixD(matrix,TMatrixD::kMultTranspose,rotmatrix));
}

string Generate_key(int ihalf, int iplane, int imodule, int iside, int idie, int isensor){
	stringstream keystream;
	keystream << ihalf << iplane << imodule << iside << idie << isensor;
	return keystream.str();
}

//#include <stdlib.h>

/**
	 * C++ version 0.4 char* style "itoa":
	 * Written by Lukás Chmela
	 * Released under GPLv3.
	 */
char* itoa(int value, char* result, int base) {
	// check that the base if valid
	char* last_char;
	if (base < 2 || base > 36) { *result = '\0'; return result; }
	char* ptr = result, *ptr1 = result, tmp_char;
	int tmp_value;
	do {
		tmp_value = value;
		value /= base;
		*ptr++ = "zyxwvutsrqponmlkjihgfedcba9876543210123456789abcdefghijklmnopqrstuvwxyz" [35 + (tmp_value - value * base)];
	} while ( value );
	// Apply negative sign
	if (tmp_value < 0) *ptr++ = '-';
	last_char = ptr;
	*ptr-- = '\0';
	//cout << last_char << endl;
	while(ptr1 < ptr) {
		tmp_char = *ptr;
		*ptr--= *ptr1;
		*ptr1++ = tmp_char;
	}
	return last_char;
}

string Generate_key_faster(int ihalf, int iplane, int imodule, int iside, int idie, int isensor){
	char key[100];
	char* ptr;
	ptr = itoa(ihalf, key, 10);
	ptr = itoa(iplane, ptr, 10);
	ptr = itoa(imodule, ptr, 10);
	ptr = itoa(iside, ptr, 10);
	ptr = itoa(idie, ptr, 10);
	ptr = itoa(isensor, ptr, 10);
    string result(key);
	//stringstream keystream;
	//keystream << ihalf << iplane << imodule << iside << idie << isensor;
	return result;
}

// simple node to construct
// a search tree for matrices with
// mother daughter relations
// up to 10 daughters can be handled
class TGeoMatrixNode {
public:
	TGeoMatrixNode () :
		fkey(-1), fmother(NULL), ftrafomatrix(NULL){
		for (unsigned char imatrix = 0; imatrix < 9; imatrix++){
			fdaughters[imatrix] = NULL;
		}
	};
	TGeoMatrixNode (signed char key, TGeoMatrixNode* mother, TGeoMatrix* trafomatrix = NULL):
		fkey(key), fmother(mother), ftrafomatrix(trafomatrix)
	{
		for (unsigned char imatrix = 0; imatrix < 9; imatrix++){
			fdaughters[imatrix] = NULL;
		}
	};
	virtual ~TGeoMatrixNode();
	void SetDaughter(signed char key,  TGeoMatrixNode* daughter);
	TGeoMatrixNode* GetDaughter (signed char key);
	void SetTrafoMatrix(TGeoMatrix* trafomatrix);
	TGeoMatrix* GetTrafoMatrix(){return ftrafomatrix;};
private:
	// key can be -1 to 9 and is the
	// position in the daughters array of the mother node
	signed char fkey;
	TGeoMatrixNode* fdaughters[10];
	TGeoMatrixNode* fmother;
	// the Geo transformation matrix into the reference system
	// of this node in respect to the mother node
	TGeoMatrix* ftrafomatrix;
};

TGeoMatrixNode::~TGeoMatrixNode (){
	delete ftrafomatrix;
	for (unsigned char imatrix = 0; imatrix < 9; imatrix++){
		delete fdaughters[imatrix];
	}
}

void TGeoMatrixNode::SetDaughter(signed char key, TGeoMatrixNode* daughter){
	if (-1 < key && key < 10){
		if (fdaughters[key] == NULL){
			fdaughters[key] = daughter;
		} else {
			cout << " Error in TGeoMatrixNode::AddDaughter: daughter for key " << key << " not empty!" << endl;
		}
	} else {
		cout << " Error in TGeoMatrixNode::AddDaughter: key " << key << " not valid!" << endl;
	}
}

TGeoMatrixNode* TGeoMatrixNode::GetDaughter (signed char key){
	if (-1 < key && key < 10){
		return fdaughters[key];
	} else {
		cout << " Error in TGeoMatrixNode::GetDaughter: key " << key << " not valid!" << endl;
		return NULL;
	}
}

void TGeoMatrixNode::SetTrafoMatrix(TGeoMatrix* trafomatrix){
	if (ftrafomatrix == NULL){
		ftrafomatrix = trafomatrix;
	} else {
		cout << " Error in TGeoMatrixNode::SetTrafoMatrix: ftrafomatrix not empty!" << endl;
	}
}

TGeoMatrixNode* froot(NULL);
TGeoMatrixNode* pGeoMatrixNode;

void AddTrafoMatrix(TGeoMatrix* trafomatrix, int ihalf = -1, int iplane = -1, int imodule = -1, int iside = -1, int idie = -1, int isensor = -1){
	if (froot == NULL){
		froot = new TGeoMatrixNode(-1, NULL, NULL);
	}
	pGeoMatrixNode = froot;
	if (ihalf == -1){
		pGeoMatrixNode->SetTrafoMatrix(trafomatrix);
	} else {
		if (pGeoMatrixNode->GetDaughter(ihalf) == NULL){
			pGeoMatrixNode->SetDaughter(ihalf, new TGeoMatrixNode(ihalf, pGeoMatrixNode, NULL));
		}
		pGeoMatrixNode = pGeoMatrixNode->GetDaughter(ihalf);
		if (iplane == -1){
			pGeoMatrixNode->SetTrafoMatrix(trafomatrix);
		} else {
			if (pGeoMatrixNode->GetDaughter(iplane) == NULL){
				pGeoMatrixNode->SetDaughter(iplane, new TGeoMatrixNode(iplane, pGeoMatrixNode, NULL));
			}
			pGeoMatrixNode = pGeoMatrixNode->GetDaughter(iplane);
			if (imodule == -1){
				pGeoMatrixNode->SetTrafoMatrix(trafomatrix);
			} else {
				if (pGeoMatrixNode->GetDaughter(imodule) == NULL){
					pGeoMatrixNode->SetDaughter(imodule, new TGeoMatrixNode(imodule, pGeoMatrixNode, NULL));
				}
				pGeoMatrixNode = pGeoMatrixNode->GetDaughter(imodule);
				if (iside == -1){
					pGeoMatrixNode->SetTrafoMatrix(trafomatrix);
				} else {
					if (pGeoMatrixNode->GetDaughter(iside) == NULL){
						pGeoMatrixNode->SetDaughter(iside, new TGeoMatrixNode(iside, pGeoMatrixNode, NULL));
					}
					pGeoMatrixNode = pGeoMatrixNode->GetDaughter(iside);
					if (idie == -1){
						pGeoMatrixNode->SetTrafoMatrix(trafomatrix);
					} else {
						if (pGeoMatrixNode->GetDaughter(idie) == NULL){
							pGeoMatrixNode->SetDaughter(idie, new TGeoMatrixNode(idie, pGeoMatrixNode, NULL));
						}
						pGeoMatrixNode = pGeoMatrixNode->GetDaughter(idie);
						if (isensor == -1){
							pGeoMatrixNode->SetTrafoMatrix(trafomatrix);
						} else {
							if (pGeoMatrixNode->GetDaughter(isensor) == NULL){
								pGeoMatrixNode->SetDaughter(isensor, new TGeoMatrixNode(isensor, pGeoMatrixNode, NULL));
							}
							pGeoMatrixNode = pGeoMatrixNode->GetDaughter(isensor);
							pGeoMatrixNode->SetTrafoMatrix(trafomatrix);
						}
					}
				}
			}
		}
	}
}

TGeoMatrix* GetTrafoMatrix(int ihalf = -1, int iplane = -1, int imodule = -1, int iside = -1, int idie = -1, int isensor = -1){
	TGeoMatrix* result(NULL);
	if (froot == NULL){
		cout << "Error: transformation matrices not loaded!" << endl;
		return NULL;
	} else {
		pGeoMatrixNode = froot;
		if (ihalf == -1){
			return pGeoMatrixNode->GetTrafoMatrix();
		} else {
			pGeoMatrixNode = pGeoMatrixNode->GetDaughter(ihalf);
			if (pGeoMatrixNode == NULL){
				cout << "Fatal: missing daughter in tree for half = " << ihalf << endl;
				return NULL;
			}
			if (iplane == -1){
				return pGeoMatrixNode->GetTrafoMatrix();
			} else {
				pGeoMatrixNode = pGeoMatrixNode->GetDaughter(iplane);
				if (pGeoMatrixNode == NULL){
					cout << "Fatal: missing daughter in tree for plane = " << iplane << endl;
					return NULL;
				}
				if (imodule == -1){
					return pGeoMatrixNode->GetTrafoMatrix();
				} else {
					pGeoMatrixNode = pGeoMatrixNode->GetDaughter(imodule);
					if (pGeoMatrixNode == NULL){
						cout << "Fatal: missing daughter in tree for module = " << imodule << endl;
						return NULL;
					}
					if (iside == -1){
						return pGeoMatrixNode->GetTrafoMatrix();
					} else {
						pGeoMatrixNode = pGeoMatrixNode->GetDaughter(iside);
						if (pGeoMatrixNode == NULL){
							cout << "Fatal: missing daughter in tree for side = " << iside << endl;
							return NULL;
						}
						if (idie == -1){
							return pGeoMatrixNode->GetTrafoMatrix();
						} else {
							pGeoMatrixNode = pGeoMatrixNode->GetDaughter(idie);
							if (pGeoMatrixNode == NULL){
								cout << "Fatal: missing daughter in tree for die = " << idie << endl;
								return NULL;
							}
							if (isensor == -1){
								return pGeoMatrixNode->GetTrafoMatrix();
							} else {
								pGeoMatrixNode = pGeoMatrixNode->GetDaughter(isensor);
								if (pGeoMatrixNode == NULL){
									cout << "Fatal: missing daughter in tree for sensor = " << isensor << endl;
									return NULL;
								}
								return pGeoMatrixNode->GetTrafoMatrix();
							}
						}
					}
				}
			}
		}
	}
	return result;
}

map<string, TGeoMatrix*> trafomatrixmap;

/*
class Tkey {
public:
	signed char half;
	signed char plane;
	signed char module;
	signed char side;
	signed char die;
	signed char sensor;
	bool operator < (const Tkey & comp) const{
		if (half < comp.half) return true;
		if (half > comp.half) return false;
		if (plane < comp.plane) return true;
		if (plane > comp.plane) return false;
		if (module < comp.module) return true;
		if (module > comp.module) return false;
		if (side < comp.side) return true;
		if (side > comp.side) return false;
		if	(die < comp.die) return true;
		if	(die > comp.die) return false;
		if	(sensor < comp.sensor) return true;
		if	(sensor >= comp.sensor) return false;
	}

	bool operator == (const Tkey & comp) const{
		return (half == comp.half) &&
			(plane == comp.plane) &&
			(module == comp.module) &&
			(side == comp.side) &&
			(die == comp.die) &&
			(sensor == comp.sensor);
	}

	Tkey (const Tkey& copy){
			half = copy.half;
			plane = copy.plane;
			module = copy.module;
			side = copy.side;
			die = copy.die;
			sensor = copy.sensor;
	}

	Tkey (int ihalf, int iplane, int imodule, int iside, int idie, int isensor){
		half = ihalf;
		plane = iplane;
		module = imodule;
		side = iside;
		die = idie;
		sensor = isensor;
	}

	Tkey (string key){
		int sign (1);
		int ikey(0); // 0,1,2,3,4,5 = ihalf, iplane, imodule, iside, idie, isensor
		int number;
		for (unsigned int ichar = 0; ichar < key.size(); ichar++){
			if (key[ichar] == '-'){
				sign = -1;
			}
			if (isdigit(key[ichar])){
				number = (key[ichar]-'0')*sign;
				sign = 1;
				if (ikey == 0) half = number;
				if (ikey == 1) plane = number;
				if (ikey == 2) module = number;
				if (ikey == 3) side = number;
				if (ikey == 4) die = number;
				if (ikey == 5) sensor = number;
				ikey++;
			}
		}
		if (ikey != 6) cout << " Error in Generate_Tkey: key string " << key << " is not valid " << endl;
	}
};*/


class Tkey {
public:
	signed char key[6];
	bool operator < (const Tkey & comp) const{
		for (unsigned char ikey = 0; ikey < 6; ikey++){
			if (key[ikey] < comp.key[ikey]) return false;
			if (key[ikey] > comp.key[ikey]) return true;
		}
		return false;
	}

	bool operator == (const Tkey & comp) const{
		for (unsigned char ikey = 0; ikey < 6; ikey++){
			if (key[ikey] != comp.key[ikey]) return false;
		}
		return true;
	}

	Tkey (const Tkey& copy){
		for (unsigned char ikey = 0; ikey < 6; ikey++){
			key[ikey] = copy.key[ikey];
		}
	}

	Tkey (int ihalf, int iplane, int imodule, int iside, int idie, int isensor){
		key[0] = ihalf;
		key[1] = iplane;
		key[2] = imodule;
		key[3] = iside;
		key[4] = idie;
		key[5] = isensor;
	}

	Tkey (string skey){
		int sign (1);
		int ikey(0); // 0,1,2,3,4,5 = ihalf, iplane, imodule, iside, idie, isensor
		int number;
		for (unsigned int ichar = 0; ichar < skey.size(); ichar++){
			if (skey[ichar] == '-'){
				sign = -1;
			}
			if (isdigit(skey[ichar])){
				number = (skey[ichar]-'0')*sign;
				sign = 1;
				key[ikey] = number;
				ikey++;
			}
		}
		if (ikey != 6) cout << " Error in Generate_Tkey: key string " << key << " is not valid " << endl;
	}
};



map<Tkey, TGeoMatrix*> trafomatrixoptmap;

#include <cstdint>
#include <omp.h>

int main(void)
{
	//TApplication myapp("myapp",0,0);
	track_finder();
	//timebased_study();
	//MickeyMouse_MC_asymetries();
	//MickeyMouse_MC_resolution();
	/*
	int nmatrices(0);
	// fill the trees
	for (int ihalf = -1; ihalf < 2; ihalf++)
		for (int iplane = -1; iplane < 0  || (ihalf > -1 && iplane < 4); iplane++)
			for (int imodule = -1; imodule < 0  || (iplane > -1 && imodule < 5); imodule++)
				for (int iside = -1; iside < 0  || (imodule > -1 && iside < 2); iside++)
					for (int idie = -1; idie < 0  || (iside > -1 && idie < 2); idie++)
						for (int isensor = -1; isensor < 0 || (idie > -1 && isensor < 3); isensor++)
			//for (unsigned int sensorID = 0; sensorID < 400 ; sensorID++)
				{
							//if (idie != -1 && isensor == -1) continue;
					string key = Generate_key_faster(ihalf, iplane, imodule, iside, idie, isensor);
					TGeoMatrix* matrix = new TGeoHMatrix(key.c_str());
					matrix->SetDx(gRandom->Uniform());
					matrix->SetDy(gRandom->Uniform());
					matrix->SetDz(gRandom->Uniform());
					matrix->RotateX(gRandom->Uniform());
					matrix->RotateY(gRandom->Uniform());
					//TGeoMatrix* matrix = new TGeoHMatrix((key + "_comp").c_str());
					AddTrafoMatrix(matrix, ihalf, iplane, imodule, iside, idie, isensor);
					trafomatrixmap[key] = matrix;
					Tkey keystruct(ihalf, iplane, imodule, iside, idie, isensor);
					trafomatrixoptmap[keystruct] = matrix;
					Tkey testkey(Generate_key_faster(ihalf, iplane, imodule, iside, idie, isensor));
					if (!(testkey == key)){
						cout << " noooo " << endl;
					}
					nmatrices++;
				}
	cout << " filled " << nmatrices << " matrices " << endl;
	cout << " number of matrices is " << trafomatrixoptmap.size() << endl;
	nmatrices = 0;
	for (map<Tkey, TGeoMatrix*>::iterator it = trafomatrixoptmap.begin(); it != trafomatrixoptmap.end(); it++){
		nmatrices++;
	}
	cout << " found " << nmatrices << endl;

	TGeoMatrix* topnode = froot->GetTrafoMatrix();

	// now the benchmark part search for the matrices in 3 different depths of the tree
	//Some benchmark tests
	vector < double > times;
	vector < string > keys;
	for (int i = 0; i < 100000; i++)
	{
		keys.clear();
		// Starting the time measurement
		double start = omp_get_wtime();
		// Computations to be measured
		//int ihalf(-1), iplane(-1), imodule(-1), iside(-1), idie(-1), isensor(-1);

		for (int ihalf = -1; ihalf < 2; ihalf++)
			for (int iplane = -1; iplane < 0  || (ihalf > -1 && iplane < 4); iplane++)
				for (int imodule = -1; imodule < 0  || (iplane > -1 && imodule < 5); imodule++)
					for (int iside = -1; iside < 0  || (imodule > -1 && iside < 2); iside++)
						for (int idie = -1; idie < 0  || (iside > -1 && idie < 2); idie++)
							for (int isensor = -1; isensor < 0 || (idie > -1 && isensor < 3); isensor++)
		//for (unsigned int sensorID = 0; sensorID < 400 ; sensorID++)
			{
			//if (idie == 1 && isensor == 0) continue;
				//string name(GetTrafoMatrix(ihalf, iplane, imodule, iside, idie, isensor)->GetName());
				//keys.push_back(name);
				TGeoMatrix* matrix1 = GetTrafoMatrix(-1, -1, -1, -1, -1, -1);
				TGeoMatrix* matrix2 = GetTrafoMatrix(ihalf, iplane, imodule, iside, -1, -1);
				TGeoMatrix* matrix3 = GetTrafoMatrix(ihalf, iplane, imodule, iside, idie, isensor);
				string name(matrix3->GetName());
				keys.push_back(name);
				//if (!matrix1 || !matrix2 || !matrix3) return TGeoHMatrix();
				TGeoHMatrix result = (((*matrix1) * (*matrix2)) * (*matrix3));
				double master[3];
				TVector3 point(1,1,1);
				point.GetXYZ(master);
				double local[3];
				result.MasterToLocal(master, local);
				//return TVector3(local);
			}

		// Measuring the elapsed time
		double end = omp_get_wtime();
		// Time calculation (in seconds)
		times.push_back(end-start);
		//cout << " time it took was " << end-start << endl;
	}

	double mean = 0;
	for (int i = 0; i < times.size(); i++){
		mean += times[i];
	}
	mean /= times.size();
	double rms = 0;
	for (int i = 0; i < times.size(); i++){
		rms += (times[i]-mean)*(times[i]-mean);
	}
	rms = sqrt(rms/(times.size()-1));

	cout << " mean to read "<< keys.size() <<" matrices from tree was " << mean << " +- " << rms << endl;

	times.clear();
	vector < string > keysf;
	for (int i = 0; i < 100000; i++)
	{
		keysf.clear();
		// Starting the time measurement
		double start = omp_get_wtime();
		// Computations to be measured
		//int ihalf(-1), iplane(-1), imodule(-1), iside(-1), idie(-1), isensor(-1);

		for (int ihalf = -1; ihalf < 2; ihalf++)
			for (int iplane = -1; iplane < 0  || (ihalf > -1 && iplane < 4); iplane++)
				for (int imodule = -1; imodule < 0  || (iplane > -1 && imodule < 5); imodule++)
					for (int iside = -1; iside < 0  || (imodule > -1 && iside < 2); iside++)
						for (int idie = -1; idie < 0  || (iside > -1 && idie < 2); idie++)
							for (int isensor = -1; isensor < 0 || (idie > -1 && isensor < 3); isensor++)
		//for (unsigned int sensorID = 0; sensorID < 400 ; sensorID++)
			{
			//if (idie == 1 && isensor == 0) continue;
				//string name(trafomatrixmap[Generate_key_faster(ihalf, iplane, imodule, iside, idie, isensor)]->GetName());
				//string name(trafomatrixoptmap[Generate_key_struct(ihalf, iplane, imodule, iside, idie, isensor)]->GetName());
				//string name(trafomatrixmap[Generate_key(ihalf, iplane, imodule, iside, idie, isensor)]->GetName());
				//string name(topnode->GetName());
				//keysf.push_back(name);
				//;
				//Tkey key2(ihalf, iplane, imodule, iside, -1, -1);
				//Tkey key3(ihalf, iplane, imodule, iside, idie, isensor);
				//TGeoMatrix* matrix1 = trafomatrixoptmap[Tkey(-1, -1, -1, -1, -1, -1)];
				//TGeoMatrix* matrix2 = trafomatrixoptmap[Tkey(ihalf, iplane, imodule, iside, -1, -1)];
				//TGeoMatrix* matrix3 = trafomatrixoptmap[Tkey(ihalf, iplane, imodule, iside, idie, isensor)];
				TGeoMatrix* matrix1 = trafomatrixmap[Generate_key(-1, -1, -1, -1, -1, -1)];
				TGeoMatrix* matrix2 = trafomatrixmap[Generate_key(ihalf, iplane, imodule, iside, -1, -1)];
				TGeoMatrix* matrix3 = trafomatrixmap[Generate_key(ihalf, iplane, imodule, iside, idie, isensor)];
				string name(matrix3->GetName());
				keysf.push_back(name);
				//if (!matrix1 || !matrix2 || !matrix3) return TGeoHMatrix();
				TGeoHMatrix result = (((*matrix1) * (*matrix2)) * (*matrix3));
				double master[3];
				TVector3 point(1,1,1);
				point.GetXYZ(master);
				double local[3];
				result.MasterToLocal(master, local);
			}

		// Measuring the elapsed time
		double end = omp_get_wtime();
		// Time calculation (in seconds)
		times.push_back(end-start);
		//cout << " time it took was " << end-start << endl;
	}

	double meanf = 0;
	for (int i = 0; i < times.size(); i++){
		meanf += times[i];
	}
	meanf /= times.size();
	double rmsf = 0;
	for (int i = 0; i < times.size(); i++){
		rmsf += (times[i]-meanf)*(times[i]-meanf);
	}
	rmsf = sqrt(rmsf/(times.size()-1));

	cout << " mean to read "<< keysf.size() <<" matrices from map was " << meanf << " +- " << rmsf << endl;

	cout << " ratio is " << meanf/mean << endl;

	// check the keys

	for (int ikey = 0; ikey < keys.size(); ikey++){
		cout << keys[ikey] << " " << keysf[ikey] << endl;
	}*/

	// mean to read 771 matrices from tree was 0.000196488 +- 2.37427e-06
	// mean to read 771 matrices from map was 0.000643847 +- 2.92558e-06
	// ratio is 3.27678 for the standard string keyed tree

	// mean to read 771 matrices from tree was 0.000153953 +- 2.40678e-06
	// mean to read 771 matrices from map was 0.000223321 +- 1.32038e-06
	// ratio is 1.45058 for an optimized struct key // bug found

	// mean to read 0 matrices from tree was 0.000132852 +- 1.43745e-06
	// mean to read 0 matrices from map was 0.000185697 +- 1.12048e-06
	// ratio is 1.39778	for an optimized struct key // bug found
	// when not storing the result of the returned value

	// mean to read 1 matrices from tree was 1.67414e-07 +- 2.96759e-08
	// mean to read 1 matrices from map was 2.85954e-07 +- 9.64492e-08
	// ratio is 1.70806 when calling only highest node

	// mean to read 1 matrices from tree was 2.26099e-07 +- 4.7714e-08
	// mean to read 1 matrices from map was 2.38775e-07 +- 5.42969e-08
	// ratio is 1.05606 when calling directly a reference to the top node

	// mean to read 0 matrices from tree was 0.000184237 +- 1.04919e-05
	// mean to read 0 matrices from map was 0.00153433 +- 3.61806e-06
	// ratio is 8.328 compared to tree with not optimized string key generator (not writing out values)

	// test with 3 matrix calculations

	// mean to read 0 matrices from tree was 0.000278152 +- 2.60515e-06
	// mean to read 0 matrices from map was 0.00408809 +- 8.61365e-06
	// ratio is 14.6973 non optimized code compared to special tree

	// mean to read 0 matrices from tree was 0.000284906 +- 2.18998e-06
	// mean to read 0 matrices from map was 0.000506343 +- 1.94691e-06
	// ratio is 1.77723 compared to optimized struct tree // bug found

	// test with 3 matrix calculations and a vector transformation

	// mean to read 0 matrices from tree was 0.00034494 +- 2.62637e-06
	// mean to read 0 matrices from map was 0.000567812 +- 2.27073e-06
	// ratio is 1.64612 compared to optimized struct tree // bug found

	// test with 3 matrix calculations and a vector transformation with non zero matrices

	// mean to read 0 matrices from tree was 0.000675651 +- 3.13462e-06
	// mean to read 0 matrices from map was 0.000905463 +- 1.92707e-06
	// ratio is 1.34013 compared to optimized struct tree // bug found

	// mean to read 0 matrices from tree was 0.000700593 +- 8.27914e-06
	// mean to read 0 matrices from map was 0.00244556 +- 4.64432e-06
	// ratio is 3.4907 compared to tree with optimized key generator

	// mean to read 0 matrices from tree was 0.000715011 +- 3.01295e-05
	// mean to read 0 matrices from map was 0.00479929 +- 5.02616e-05
	// ratio is 6.71219 compared to the non optimized tree

	// mean to read 0 matrices from tree was 0.000702061 +- 6.80643e-06
	// mean to read 0 matrices from map was 0.00145913 +- 1.83538e-05
	// ratio is 2.07835 compared to optimized class tree

	// mean to read 0 matrices from tree was 0.000751482 +- 4.51953e-06
	// mean to read 0 matrices from map was 0.00471662 +- 1.26561e-05
	// ratio is 6.27642 compared to non optimized tree

	// mean to read 0 matrices from tree was 0.000701484 +- 3.54751e-06
	// mean to read 0 matrices from map was 0.00249213 +- 4.55697e-06
	// ratio is 3.55265 compared to tree with optimized key generator

	// ************ char array as key ***********
	// mean to read 771 matrices from tree was 0.00106938 +- 4.93655e-06
	// mean to read 771 matrices from map was 0.00262172 +- 1.54522e-05
	// ratio is 2.45162 old map with faster key generator

	// mean to read 771 matrices from tree was 0.00100945 +- 5.48069e-06
	// mean to read 771 matrices from map was 0.00183734 +- 9.23518e-06
	// ratio is 1.82014 for map with class as key (char[6])

	// mean to read 771 matrices from tree was 0.00102953 +- 1.39849e-05
	// mean to read 771 matrices from map was 0.00191463 +- 3.80041e-06
	// ratio is 1.85972 same as above but defined as struct instead of class

	// mean to read 771 matrices from tree was 0.00092196 +- 9.78848e-06
	// mean to read 771 matrices from map was 0.00165397 +- 1.25088e-05
	// ratio is 1.79397 for map with class as key and differentiated members

	// mean to read 771 matrices from tree was 0.000995503 +- 4.97537e-06
	// mean to read 771 matrices from map was 0.00502832 +- 1.9187e-05
	// ratio is 5.05103 for no optimization



	//cout << " tree is deleted " << endl;
	/*
	//Some benchmark tests
	vector < double > times;
	vector < string > keys;
	for (int i = 0; i < 1000; i++)
	{
		keys.clear();
		// Starting the time measurement
		double start = omp_get_wtime();
		// Computations to be measured

		for (int ihalf = -1; ihalf < 2; ihalf++)
			for (int iplane = -1; iplane < 4; iplane++)
				for (int imodule = -1; imodule < 5; imodule++)
					for (int iside = -1; iside < 2; iside++)
						for (int idie = -1; idie < 2; idie++)
							for (int isensor = -1; isensor < 3; isensor++)
		//for (unsigned int sensorID = 0; sensorID < 400 ; sensorID++)
			{
			if (idie == 1 && isensor == 0) continue;
				keys.push_back(Generate_key(ihalf, iplane, imodule, iside, idie, isensor));
			}

		// Measuring the elapsed time
		double end = omp_get_wtime();
		// Time calculation (in seconds)
		times.push_back(end-start);
		//cout << " time it took was " << end-start << endl;
	}

	double mean = 0;
	for (int i = 0; i < times.size(); i++){
		mean += times[i];
	}
	mean /= times.size();
	double rms = 0;
	for (int i = 0; i < times.size(); i++){
		rms += (times[i]-mean)*(times[i]-mean);
	}
	rms = sqrt(rms/(times.size()-1));

	cout << " mean to generate "<< keys.size() <<" keys was " << mean << " +- " << rms << endl;

	times.clear();
	vector < string > keysf;
	for (int i = 0; i < 1000; i++)
	{
		keysf.clear();
		// Starting the time measurement
		double start = omp_get_wtime();
		// Computations to be measured

		for (int ihalf = -1; ihalf < 2; ihalf++)
			for (int iplane = -1; iplane < 4; iplane++)
				for (int imodule = -1; imodule < 5; imodule++)
					for (int iside = -1; iside < 2; iside++)
						for (int idie = -1; idie < 2; idie++)
							for (int isensor = -1; isensor < 3; isensor++)
		//for (unsigned int sensorID = 0; sensorID < 400 ; sensorID++)
			{
			if (idie == 1 && isensor == 0) continue;
				keysf.push_back(Generate_key_faster(ihalf, iplane, imodule, iside, idie, isensor));
			}

		// Measuring the elapsed time
		double end = omp_get_wtime();
		// Time calculation (in seconds)
		times.push_back(end-start);
		//cout << " time it took was " << end-start << endl;
	}

	mean = 0;
	for (int i = 0; i < times.size(); i++){
		mean += times[i];
	}
	mean /= times.size();
	rms = 0;
	for (int i = 0; i < times.size(); i++){
		rms += (times[i]-mean)*(times[i]-mean);
	}
	rms = sqrt(rms/(times.size()-1));

	cout << " mean to generate "<< keysf.size() <<" keys faster was " << mean << " +- " << rms << endl;

	// check the keys
	for (int ikey = 0; ikey < keys.size(); ikey++){
		cout << keys[ikey] << " " << keysf[ikey] << endl;

	}

	// not optimized result:
	// mean to generate 2970 keys was 0.00314231 +- 2.03523e-05
	// without the << operator part in the string twice as fast:
	//  mean to generate 2970 keys was 0.00164725 +- 9.41036e-06
	// a direct conversion via sequences of itoa
	// mean to generate 2970 keys was 0.000637954 +- 9.02265e-06
*/

	/*
	//MickeyMouse_MC_resolution();

	unsigned char chopped_number [4];
	for (int i = 0; i < 10; i++){
		// generate random 8 bit number
		chopped_number[0]=rand()%256;
		chopped_number[1]=rand()%256;
		chopped_number[2]=rand()%256;
		chopped_number[3]=rand()%256;
		cout << (int)chopped_number[0] << " " << (int)chopped_number[1] << " " << (int)chopped_number[2] << " " << (int)chopped_number[3] << endl;
		cout << ((int)chopped_number[0]*24)+((int)chopped_number[2]*16)+((int)chopped_number[1]*8)+((int)chopped_number[3]*0) << endl;
		uint32_t number = (chopped_number[0]<<24) | (chopped_number[1]<<16) | (chopped_number[2]<<8 | (chopped_number[0]<<0));
		cout << number << endl << endl;
	}*/

	/*
	unsigned char tmp [4];
	unsigned char *dataaddr = &Data[ReadPointer];
	tmp[0]  = *(dataaddr);
	tmp[1]  = *(dataaddr+1);
	tmp[2]  = *(dataaddr+2);
	tmp[3]  = *(dataaddr+3);
	ReadPointer = (ReadPointer + 4)%BufferSize;
	return (((uint32_t)tmp[0] << 24) | ((uint32_t) tmp[1] << 16) | ((uint32_t) tmp[2] << 8)  | ((uint32_t) tmp[3] << 0));
*/
	/*
	string compstr("lmd_testbox_000");
	cout << compstr.size() << endl;
	const char* anotherstr = "lmd";
	cout << compstr.compare(0, 3, anotherstr) << endl;
	*/
	/*
	cout << "" << endl; // prints
	float testfloat = 1242403838;
	double testdouble = 1242403838;//(double) testfloat;
	cout << testfloat << endl;
	cout << testdouble << endl;
	int testint = floor(testdouble);
	cout << (int) testfloat  << endl;
	cout << (int) testdouble << endl;
	cout << testint << endl;
	return 0;
*/
	/*
	TMatrixD matrix(Get_matrix_for_rotation());
	TMatrixD& rotated_matrix1 = Transform_global_to_lmd_local_ref(matrix);
	const TMatrixD& rotated_matrix2 = Transform_global_to_lmd_local_val(matrix);
	const TMatrixD& rotated_matrix3 = Transform_global_to_lmd_local_ana(matrix);

	rotated_matrix1.Print();
	rotated_matrix2.Print();
	rotated_matrix3.Print();
*/
	//TApplication myapp("myapp",0,0);
	//box_fourier_transform();
	//Construct_plane();
	//MickeyMouse_MC_resolution();
	/*
	{
		cout << " test 1" << endl;
		const testclass& myclass = get_testclass();
		//myclass.~testclass();
		cout << " got the object " << endl;
		//delete myclass;
		//myclass.~testclass();
	}
	cout << " done " << endl;
	{
		cout << " test 2" << endl;
		testclass myclass = get_testclass();
		cout << " got the object " << endl;
	}
	cout << " done " << endl;
	{
		cout << " test 3" << endl;
		testclass myclass(get_testclass());
		cout << " got the object " << endl;
	}
	cout << " done " << endl;*/


	//myapp.Run();
	return 0;
}


/*
double T = 0.5; // mus 150 m buch length over 3*10^8 m/s speed of light
double L = 1.91; // mus 574 m HESR length over 3*10^8 m/s speed of light

Double_t periodic_box_function(Double_t *x, Double_t *par){
	double xx = x[0]+T/2.; // center the box around zero
	xx = xx - L * floor(xx/L);
	if (0. < xx && xx < T) {return 1.;}
	else {return 0.;}
}

//Double_t fourier_harmonic(Double_t *x, Double_t *par){
//	double result = 0.;
//	double xx = x[0]+T/2.; // center the periodic signal around zero
//	result = ()
//	return result;
//}

void box_fourier_transform(){
	TCanvas* canvas = new TCanvas("canvas", "canvas", 800, 600);
	const int nbins = 1000;
	TF1* func_box = new TF1("func_box", periodic_box_function, -L/2.*3., +L/2.*3.);
	TH1F* hist_box = new TH1F("hist_box","hist_box", nbins, -2*L, 2*L);
	for (int ibin = 1; ibin <= nbins; ibin++){
		hist_box->SetBinContent(ibin, func_box->Eval(hist_box->GetBinCenter(ibin)));
	}
//
//	func_box->SetNpx(1000);
	hist_box->Draw("");
	hist_box->GetXaxis()->SetTitle("time / #mus");
	hist_box->GetYaxis()->SetTitle("amplitude");
	canvas->Print("box_fourier_transform.pdf");

	TH1 *hm = NULL;
	TVirtualFFT::SetTransform(0);
	hm = hist_box->FFT(hm, "MAG");
	if (hm){
		hm->SetTitle("Magnitude of the 1st transform");
		hm->Draw();

		canvas->Print("box_fourier_transform_mag.pdf");
	}

}

// to run type
// root -l Construct_plane.C+ -q

#include "TCanvas.h"
#include "TStyle.h"
#include "TROOT.h"
#include "TLine.h"
#include "TPolyLine.h"
#include "TMath.h"
#include "TH2.h"
#include "TGraph.h"
#include <iostream>
#include "TLatex.h"
#include <sstream>
#include <fstream>

using namespace std;

// two circles to fit electronics in between
float outer_radius = 90; // mm
float inner_radius = 36; // mm

// current design from report:
// ------------------------- dead_height (electronics)
// -------------------------
// |       |       |       |
// |       |       |       | height (pixel area)
// |       |       |       |
// -------------------------
//   width   width   width
//
const float detector_width  = 20; // mm
const float detector_height = 20; // mm
const float detector_dead_height = 0.5; // mm
const float detector_dead_width  = 0.; // mm
const float detector_thickness = 0.005; // mm (not used)
const int n_detectors_in_row = 3; // number of detectors per silicon waver
const int n_detectors_in_col = 1; // currently only one
const float detectors_spacing = 1.; // spacing between silicon wavers (guess)
const float pixel_size = 0.05;
const int npixels_per_detector = floor(detector_width/pixel_size * detector_height/pixel_size);
const float power_consumption = 7; // mW / mm˛

const float pi = 3.141;

// TArc made me problems drawing it not opaque with Fill option 4000
// therefore I've written my own method to create a circle
TPolyLine* Get_PolyLine_circle(float x,float y, float r);

// Draw the circles according to the inner and outer radius as given
// A coordinate system is drawn as well
// A reference to the graph that is drawn is given
TGraph* Draw_Circles();

// both arrays of length 5 contain the edges of the modules
// Those coordinates are rotated around the origin by the given angle
void Rotate_module(double* xpoints, double* ypoints, double angle);

void Construct_plane(){
	gROOT->SetStyle("Plain");
	// plotting if requested
	TCanvas* canvas_visualization = new TCanvas("visualization", "visualization", 1000, 1000);
	//gPad->Range(-outer_radius*(1+0.1), -outer_radius*(1+0.1), outer_radius*(1+0.1), outer_radius*(1+0.1));
	TGraph* world_coord_system = Draw_Circles();
	// write some values out
	ofstream outputfile("module_positions.dat");
	// place the parts starting from the mean y and x line as the detector is supposed to be divided into 4 parts
	// start from the  edge
	double module_width  = (detector_width + detector_dead_width)*n_detectors_in_row;
	double module_height = (detector_height + detector_dead_height)*n_detectors_in_col;
	outputfile << " module size " << module_width << " x " << module_height << " mm^2 " << endl;
	float y = 0.;
	int n_pixels = 0;
	int n_modules = 0;
	while (y < outer_radius*(1+0.1)) {
		float x = 0;
		// place it next to the hole
		if (y < inner_radius){
			x = sqrt(inner_radius*inner_radius - y*y);
		}
		while (x < outer_radius*(1+0.1)) {
			// create the 4 square coordinates of the silicon waver	in the first quater of the coordinate system
			double xpoints[5] = {x, x+module_width, x+module_width, x, x};
			double ypoints[5] = {y, y, y+module_height, y+module_height, y};
			// check the first sensor to be within the two circles
			double xmean = x+detector_width/2.;
			double ymean = y+detector_height/2.;
			double r2mean = xmean*xmean + ymean*ymean;
			if ((inner_radius*inner_radius <= r2mean && r2mean <= outer_radius*outer_radius)){
				// place always 4 parts in each quarter
				TPolyLine *pline_module = new TPolyLine(5,xpoints, ypoints);
				pline_module->Draw();
				outputfile << (xpoints[1]+xpoints[0])/2. << ", " << (ypoints[2]+ypoints[1])/2. << endl;
				// mirror on y
				for (int i = 0; i < 5; i++) ypoints[i] = -ypoints[i];
				pline_module = new TPolyLine(5,xpoints, ypoints);
				pline_module->Draw();
				outputfile << (xpoints[1]+xpoints[0])/2. << ", " << (ypoints[2]+ypoints[1])/2. << endl;
				// mirror on x
				for (int i = 0; i < 5; i++) xpoints[i] = -xpoints[i];
				pline_module = new TPolyLine(5,xpoints, ypoints);
				pline_module->Draw();
				outputfile << (xpoints[1]+xpoints[0])/2. << ", " << (ypoints[2]+ypoints[1])/2. << endl;
				// mirror on y
				for (int i = 0; i < 5; i++) ypoints[i] = -ypoints[i];
				pline_module = new TPolyLine(5,xpoints, ypoints);
				pline_module->Draw();
				outputfile << (xpoints[1]+xpoints[0])/2. << ", " << (ypoints[2]+ypoints[1])/2. << endl;
				//gPad->Update();
				n_modules+=4;
				n_pixels += npixels_per_detector*n_detectors_in_row*n_detectors_in_col*4;
			}
			// move to next row
			x += module_width+detectors_spacing;
		}
		// move to next col
		y += module_height+detectors_spacing;
	}
	outputfile.close();
	stringstream sline;
	TLatex text;
	text.SetTextSize(world_coord_system->GetXaxis()->GetLabelSize()*0.4);
	text.SetTextAlign(22);
	sline << "spacing = " << detectors_spacing << " mm";
	text.DrawText(0, 15, sline.str().c_str());
	sline.str("");
	sline << "mod width = " << detector_width << " + " << detector_dead_width << " mm X " << n_detectors_in_row;
	text.DrawText(0, 10, sline.str().c_str());
	sline.str("");
	sline << "mod height = " << detector_height << " + " << detector_dead_height << " mm X " << n_detectors_in_col;
	text.DrawText(0, 5, sline.str().c_str());
	sline.str("");
	sline << "in total " << n_pixels << " pixel";
	text.DrawText(0, 0, sline.str().c_str());
	sline.str("");
	sline << "in " << n_modules << " modules";
	text.DrawText(0, -5, sline.str().c_str());
	sline.str("");
	double total_power_consumption = n_modules*module_width*module_height*power_consumption/1000.;
	sline << "est. power consump. = ";
	text.DrawText(0, -10, sline.str().c_str());
	sline.str("");
	sline << total_power_consumption << " W";
	text.DrawText(0, -15, sline.str().c_str());
	sline.str("");
	canvas_visualization->Print("visualization_1.pdf");

// try second arrangement in a circle
	canvas_visualization->Clear();
	delete world_coord_system;
	world_coord_system = Draw_Circles();
	// calculate how many silicon modules are needed to cover the radial distance
	// construct such an array
	double phi = 0.;
	double delta_phi = 0.00001;
	while ((- pi/2. + delta_phi) < pi/2.){
		phi = - pi/2. + delta_phi;
		delta_phi = 0;
		double r = inner_radius;
		while (r < outer_radius){
		double x = r;
		double y = 0.;
		double xpoints[5] = {x, x+module_width, x+module_width, x, x};
		double ypoints[5] = {y, y, y+module_height, y+module_height, y};
		Rotate_module(xpoints, ypoints, phi);
		if (delta_phi == 0){
			// the next module has to be rotated as far as it does not collide
			// with the lower left edge of the previous module
			// create a normalized vector from origin to that point
			// vect_x is not needed since x component of the reference vector is 0
			//double vect_x = xpoints[3]/(sqrt(xpoints[3]*xpoints[3]+ypoints[3]*ypoints[3]));
			double vect_norm = (sqrt(xpoints[3]*xpoints[3]+ypoints[3]*ypoints[3]));
			double vect_y = ypoints[3]/vect_norm;
			// calculate the angle between pi/2. and that vector
			if (xpoints[3]<0.) delta_phi = pi;
			else delta_phi = acos(-vect_y) + detectors_spacing/vect_norm;
			// the latter expression is a replacement of tan(spacing/norm) as the angle is small
		}
		if ((- pi/2. + delta_phi) < pi/2.){//(inner_radius*inner_radius <= r2mean && r2mean <= outer_radius*outer_radius)){
			// place always 4 parts in each quarter
			TPolyLine *pline_module = new TPolyLine(5,xpoints, ypoints);
			pline_module->Draw();
			//outputfile << (xpoints[1]+xpoints[0])/2. << ", " << (ypoints[2]+ypoints[1])/2. << endl;
			// mirror on y
			for (int i = 0; i < 5; i++) ypoints[i] = -ypoints[i];
			//pline_module = new TPolyLine(5,xpoints, ypoints);
			//pline_module->Draw();
			//outputfile << (xpoints[1]+xpoints[0])/2. << ", " << (ypoints[2]+ypoints[1])/2. << endl;
			// mirror on x
			for (int i = 0; i < 5; i++) xpoints[i] = -xpoints[i];
			pline_module = new TPolyLine(5,xpoints, ypoints);
			pline_module->Draw();
			//outputfile << (xpoints[1]+xpoints[0])/2. << ", " << (ypoints[2]+ypoints[1])/2. << endl;
			// mirror on y
			for (int i = 0; i < 5; i++) ypoints[i] = -ypoints[i];
			//pline_module = new TPolyLine(5,xpoints, ypoints);
			//pline_module->Draw();
			//outputfile << (xpoints[1]+xpoints[0])/2. << ", " << (ypoints[2]+ypoints[1])/2. << endl;
			//gPad->Update();
			n_modules+=4;
			n_pixels += npixels_per_detector*n_detectors_in_row*n_detectors_in_col*4;
		}
		gPad->Update();
		sleep(1);
		r += module_width + detectors_spacing;
		}
	}

	canvas_visualization->Print("visualization_2.pdf");
}

TGraph* Draw_Circles(){
	double xpoints[4] = {-outer_radius*(1+0.1), -outer_radius*(1+0.1),  outer_radius*(1+0.1), outer_radius*(1+0.1)};
	double ypoints[4] = {-outer_radius*(1+0.1),  outer_radius*(1+0.1), -outer_radius*(1+0.1), outer_radius*(1+0.1)};
	TGraph* world_coord_system = new TGraph(4, xpoints, ypoints);
	world_coord_system->SetNameTitle("world_coord_system", "world coordinate system");
	world_coord_system->Draw("AP");
	world_coord_system->GetXaxis()->SetTitle("X/mm");
	world_coord_system->GetYaxis()->SetTitle("Y/mm");
	TPolyLine *arc_outer = Get_PolyLine_circle(0,0,outer_radius);
    	arc_outer->SetLineColor(1);
   	 arc_outer->SetLineWidth(2);
	arc_outer->Draw();
	TPolyLine *arc_inner = Get_PolyLine_circle(0,0,inner_radius);
   	 arc_inner->SetLineColor(1);
    	arc_inner->SetLineWidth(2);
	arc_inner->Draw();
	gPad->Update();
	return world_coord_system;
}

TPolyLine* Get_PolyLine_circle(float x,float y, float r){
    const int nsegments = 360;
	Double_t _x[nsegments+1] = {0.,}; // +1 : Endpoint is Startpoint
    Double_t _y[nsegments+1] = {0.,};
    for (int i = 0; i <= nsegments; i++){
    	double angle = 2*pi*((double) i)/((double) nsegments);
    	_x[i] = r*sin(angle)+x;
    	_y[i] = r*cos(angle)+y;
    }
    TPolyLine *pline_circle = new TPolyLine(nsegments+1,_x,_y);
    return pline_circle;
}

void Rotate_module(double* xpoints, double* ypoints, double angle){

	for (int i = 0; i < 5; i++){
		double x = xpoints[i];
		double y = ypoints[i];
		xpoints[i] = x * cos(angle) - y * sin(angle);
		ypoints[i] = x * sin(angle) + y * cos(angle);
	}
}*/