//
//- Emiliano Mocchiutti
//
//  Version 1.01 (2005-11-18)
//
//
//   Changelog:
//
//   1.00 - 1.01 (2005-11-18): Use pathtocalibration routine.
//
//   0.00 - 1.00 (2005-09-09): Working.
//
//   0.00 (2005-07-25): Created.
//
#include <fstream>
#include <TTree.h>
#include <TClassEdit.h>
#include <TObject.h>
#include <TList.h>
#include <TSystem.h>
#include <TSystemDirectory.h>
#include <TString.h>
#include <TFile.h>
#include <TClass.h>
#include <TCanvas.h>
#include <TH1.h>
#include <TH1F.h>
#include <TH2D.h>
#include <TLatex.h>
#include <TPad.h>
#include <TGraph.h>
#include <TF1.h>
//
#include <string>
#include <ccal2struct.h>
#include <ctrkstruct.h>
#include <ctrkinclude.h>
//
#if !defined (__CINT__)				
extern  const char *pathtocalibration();
extern  void creadB(const char []);
extern  void ctrack(int , Double_t [], Double_t [], Double_t [], Double_t [], int &);
#include <event/PamelaRun.h>
#include <event/physics/calorimeter/CalorimeterEvent.h>
#include <event/physics/trigger/TriggerEvent.h>
#include <event/CalibCalPedEvent.h>
#include <caloclasses.h>
#endif
//
#include <CaloFunctions.h>
//

void CaloTrackerAlignmentRT(TString filename){
    gDirectory->GetList()->Delete();
#if defined (__CINT__)	
    emicheckLib();
    stringstream llib;
    const char *pamlib = gSystem->ExpandPathName("$PAM_LIB");
    llib << pamlib << "/caloclasses_h.so"; 
    gSystem->Load(llib.str().c_str());
    stringstream lolib;
    lolib.str("");
    lolib << pamlib << "/UTpathtoc_c.so";
    gSystem->Load(lolib.str().c_str());
#endif
    const char* startingdir = gSystem->WorkingDirectory();
    TString path;
    stringcopy(path,startingdir);
    //
    //
    Long64_t nevents = 0;
    Long64_t trnevents = 0;
    //
    printf("\n\n TRACKER: loading shared libraries...\n");
#if defined (__CINT__)				
    char *pam_lib=gSystem->Getenv("PAM_LIB");
    stringstream libload;
    libload.str("");
    libload << pam_lib << "/libtrack";
    gSystem->Load(libload.str().c_str());
    libload.str("");
    libload << pam_lib << "/libtrack_C";
    gSystem->Load(libload.str().c_str());
    libload.str("");
    libload << pam_lib << "/libreadb2maps";
    gSystem->Load(libload.str().c_str());
    libload.str("");
    libload << pam_lib << "/libreadb2maps_C";
    gSystem->Load(libload.str().c_str());
#endif
    printf(" ...done!\n\n");
    //
    // Read the magnetic field map (two files), here it goes the path AND name of one of them.
    //
    printf("\n\n TRACKER: loading the magnetic field maps...\n\n\n");			 
    const char *pam_calib = pathtocalibration();
    //   const char *pam_calib=gSystem->Getenv("PAM_CALIB");			
    stringstream bdir;
    bdir.str("");
    bdir << pam_calib << ".";
    creadB((char *)bdir.str().c_str());
    //
    printf(" ...done! \n");
    //    
    Float_t rig;
    //
    Float_t sigmacalo = 10.;
    //
    Float_t piano[22];
    piano[0] = 0.;
    for ( Int_t ee = 1; ee<22; ee++){
	if ( ee%2 ){
	    piano[ee] = piano[ee-1] - 8.42;
	} else {
	    piano[ee] = piano[ee-1] - 9.76;
	};
    };
    //
    Float_t zoff = 0.; // 260
    //
    Float_t chi2 = 0.;
    Float_t chi2nev = 0.;
    //
    Float_t xalign0 = -120.;
    Float_t yalign0 = -120.;
    Float_t zalign0 = -260.; // 0
    Float_t phi0 = 0.;
    Float_t the0 = 0.;
    Float_t psi0 = 0.;
    Float_t dxalig = 100.; //50
    Float_t dyalig = 100.; //50
    Float_t dzalig = 140.; //80
    //
    Float_t dpsi = 0.20;
    Float_t dthe = 0.20;
    Float_t dphi = 0.20;
    //
    TString lfile;
    TString sfile;
    TString file2f;
    TString pfile;
    TString fififile;
    //
    Float_t minxalign;
    Float_t minyalign;
    Float_t minzalign;
    Float_t minpsialign;
    Float_t minthealign;
    Float_t minphialign;
    Float_t minxalign0;
    Float_t minyalign0;
    Float_t minzalign0;
    Float_t minpsi0;
    Float_t minthe0;
    Float_t minphi0;
    Int_t minnxa;
    Int_t minnza;
    Int_t minnya;
    Int_t minphi;
    Int_t minpsi;
    Int_t minthe;
    //
    const int npia = 0;
    TFile *hfile = 0;
    //
    //
    const string myfil = (const char *)filename;
    Int_t myposiz = myfil.find("dw_");    
    if ( myposiz == -1 ) {
      myposiz = myfil.find("DW_");
    };
    TString fileid;
    stringcopy(fileid,filename,myposiz+3,myposiz+13);
    //
    //
    stringstream namefile;
    TTree *tree;
    const char *pattho = path;
    namefile.str("");
    namefile << pattho << "/CaloTrkAlignData_";
    namefile << fileid.Data()  << ".root";
    printf(" Save file %s\n",namefile.str().c_str());
    hfile = new TFile(namefile.str().c_str(),"RECREATE","Calorimeter-Tracker alignment data");
    tree = new TTree("CaloTrkAlign","Calorimeter-Tracker alignment data");
    tree->Branch("chi2",&chi2,"chi2/F");
    tree->Branch("chi2nev",&chi2nev,"chi2nev/F");
    tree->Branch("minnxa",&minnxa,"minnxa/I");
    tree->Branch("minnya",&minnya,"minnya/I");
    tree->Branch("minnza",&minnza,"minnza/I");
    tree->Branch("minpsi",&minpsi,"minpsi/I");
    tree->Branch("minphi",&minphi,"minphi/I");
    tree->Branch("minthe",&minthe,"minthe/I");
    tree->Branch("minxalign",&minxalign,"minxalign/F");
    tree->Branch("minyalign",&minyalign,"minyalign/F");
    tree->Branch("minzalign",&minzalign,"minzalign/F");
    tree->Branch("minpsialign",&minpsialign,"minpsialign/F");
    tree->Branch("minphialign",&minphialign,"minphialign/F");
    tree->Branch("minthealign",&minthealign,"minthealign/F");
    //
    const int ngrid = 11;
    Float_t togl = ((float)ngrid-1.)/2.;
    //
    struct Tracklev2 trk;
    struct CaLevel2 calo;
    //
    Int_t i = -1;    
    Int_t itr = -1;
    Int_t syncro = 1;
    Int_t pktnum = 0;
    Int_t obt = 0;    
    //
    // open files
    //    
    TFile *caloFile = 0;
    TFile *trackerFile = 0;
    TTree *tr = 0;
    TTree *otr = 0;
    //
    caloFile = emigetFile(filename, "Physics.Level2", "Calorimeter");
    if ( !caloFile ){
    	printf("\n ERROR: no calorimeter file! \n\n");
    	return;
    };
    trackerFile = emigetFile(filename, "Physics.Level2", "Tracker");
    if ( !trackerFile ){
	printf("\n ERROR: no tracker file! \n\n");
	caloFile->Close();
	return;
    }
    //
    //Takes the tree of the header file
    //    
    tr = (TTree*)trackerFile->Get("TrkLevel2");
    settrklev2(tr,trk);
    //
    otr = (TTree*)caloFile->Get("CaloLevel2");
    setcalolevel2(otr,calo);
    //
    trnevents = tr->GetEntries();
    nevents = otr->GetEntries();
    const int lung = nevents;
    Int_t goodev[lung];
    Float_t mxtb[11][lung];
    Float_t mytb[11][lung];
    //
    i = -1;    
    itr = -1;
    syncro = 1;
    pktnum = 0;
    obt = 0;
    printf(" Initializing... \n");
    Int_t npoint = 44;
    Int_t goodevno = 0;
    while ( i < (nevents-1) ){    
	//
	if ( i%1000 == 0 && i > 0 ) printf(" %iK \n",i/1000);	
	//
	// look for tracker data
	//
	itr++;
	i++;
	syncro = 1;
	//
    trkcalosync0: 
	//
	// check if we have tracker data
	//
	if ( i >= nevents ){
	    printf(" WARNING: no more calorimeter data.\n");
	    goto skipthisev0;
	};
	if ( itr > trnevents ){
	    printf(" WARNING: no more tracker data.\n");
	    goto skipthisev0;
	};
	//
	// retrieve tracker informations
	//
	tr->GetEntry(itr);
	//
	// check synchronization tracker and calorimeter informations:
	//
	otr->GetEntry(i);
	//
	pktnum = calo.pkt_num;
	obt = calo.OBT;
	if ( pktnum != trk.pkt_num || obt != trk.obt ){
	    if ( pktnum > trk.pkt_num || obt > trk.obt ){
		itr++;
		syncro = 0;		
		goto trkcalosync0;
	    };
	    if ( pktnum < trk.pkt_num || obt < trk.obt ){
		i++;
		syncro = 0;
		goto trkcalosync0;
	    };
	};	
	//
	// here we have synchronized data
	//
	if ( !syncro ) {
	    syncro = 1;
	};       
	//
	// here we go!
	//
	goodev[i] = 0;
	//
	// select events with one track...
	//
	if ( trk.ntrk == 1 && trk.good2 ){
	    //
	    // a good chi2 from tracker and rigidity greater than 1 GV (notice correction in deflection)
	    //
	    rig = 1./(trk.al[0][4]);
	    if ( trk.chi2[0] < 25 &&  (rig > 0.5 || rig < -0.5) ){
		//
		// a good fit of the track from calorimeter information
		//
		if ( calo.npcfit[0]>15 && calo.npcfit[1]>15 && calo.varcfit[0]<20. && calo.varcfit[1]<20. && calo.good == 1 ){ //era 19!
		    //
		    // determine the position of the track according to tracker:
		    //
		    Double_t zin[22];
		    Double_t xtb[22];
		    Double_t ytb[22];
		    Double_t al_p[5];
		    Int_t ifail;
		    Int_t ee;
		    Int_t eb;
		    //			
		    minzalign0 = zalign0;
		    minzalign = 0.;
		    //
		    ee = 0;
		    eb = 11;
		    //
		    for ( Int_t ea = 0; ea < 11; ea++) {
			//
			eb--;
			minzalign = minzalign0 + ((float)eb - togl)*(dzalig/((float)ngrid+1.));
			//
			if ( ea < 5 ) al_p[ea] = (Double_t)trk.al[0][ea];
		    	xtb[ee] = 0.;
		    	ytb[ee] = 0.;					    
			zin[ee] = (piano[npia] + minzalign)/10.;
		    	ee++;
		    	xtb[ee] = 0.;
		    	ytb[ee] = 0.;
		    	zin[ee] = (piano[npia] - 5.1 + minzalign)/10.;
		    	ee++;
		    };
		    //
		    npoint = 22;
		    ifail = 0;
		    //
		    ctrack(npoint,zin,xtb,ytb,al_p,ifail);	
		    //			
		    // here alignment factors for x and y
		    //
		    if ( !ifail ){
			goodev[i] = 1;
			goodevno++;
			//				    
			minzalign0 = zalign0;
			minzalign = 0.;
			ee = 21;
			for ( Int_t ea = 0; ea < 11 ; ea++) {			    
			    mxtb[ea][i] = (float)xtb[ee] * 10.;					
			    ee--;
			    mytb[ea][i] = (float)ytb[ee] * 10.;
			    ee--;
			};
			//
		    };
		};
	    };
	};
	//
    skipthisev0:
	continue;
    };    
    printf(" ...done! \n");
    printf(" %i good events \n",goodevno);
    //
    Int_t ee = 0;
    Float_t the;
    Float_t psi;
    Float_t phi;   
    //
    Float_t xcbar;
    Float_t ycbar;
    Float_t xtbar;
    Float_t ytbar;
    //
    //    for (Int_t it = 0; it < 3; it++){
    for (Int_t it = 0; it < 1; it++){
	//	printf(" iteraction %i \n",it);
	//
	if ( it == 0 ){
	    minxalign0 = xalign0;
	    minyalign0 = yalign0;
	    minzalign0 = zalign0;
	    minpsi0 = psi0;
	    minthe0 = the0;
	    minphi0 = phi0;
	} else {
	    minxalign0 = minxalign;
	    minyalign0 = minyalign;
	    minzalign0 = minzalign;
	    minpsi0 = minpsialign;
	    minthe0 = minthealign;
	    minphi0 = minphialign;
	};
	//
	for (Int_t nthe = 0; nthe<ngrid; nthe++){
	    //for (Int_t nthe = 5; nthe<6; nthe++){
	    //for (Int_t nthe = 4; nthe<7; nthe++){
	    minthealign = minthe0 + ((float)nthe - togl)*(dthe/((float)ngrid+1.));
	    //
	    for (Int_t nphi = 0; nphi<ngrid; nphi++){
	    //for (Int_t nphi = 4; nphi<5; nphi++){
	    //for (Int_t nphi = 5; nphi<6; nphi++){
		//for (Int_t nphi = 4; nphi<7; nphi++){
		minphialign = minphi0 + ((float)nphi - togl)*(dphi/((float)ngrid+1.));
		//
		for (Int_t npsi = 0; npsi<ngrid; npsi++){
		    //for (Int_t npsi = 5; npsi<6; npsi++){
		    //for (Int_t npsi = 4; npsi<7; npsi++){
		    minpsialign = minpsi0 + ((float)npsi - togl)*(dpsi/((float)ngrid+1.));
		    //
		    for (Int_t nxa = 0; nxa<ngrid; nxa++){
		    //for (Int_t nxa = 4; nxa<5; nxa++){
			//	    printf(" nxa %i \n",nxa);
			minxalign = minxalign0 + ((float)nxa - togl)*(dxalig/((float)ngrid+1.));
			//printf("minxalign[%i][%i] = %f\n",g,h,minxalign[g][h]);
			//for (Int_t nya = 6; nya<7; nya++){
			for (Int_t nya = 0; nya<ngrid; nya++){
			    //		printf(" nya %i \n",nya);
			    minyalign = minyalign0 + ((float)nya - togl)*(dyalig/((float)ngrid+1.));
			    //printf("minyalign[%i][%i] = %f\n",g,h,minyalign[g][h]);
			    //for (Int_t nza = 5; nza<6; nza++){
			    //for (Int_t nza = ngrid ; nza>=0; nza--){
			    for (Int_t nza = 0 ; nza < ngrid; nza++){
				//
				minzalign = minzalign0 + ((float)nza - togl)*(dzalig/((float)ngrid+1.));
				printf(" Step %i xa %f ya %f za %f phi %f the %f psi %f \n",it,minxalign,minyalign,minzalign,minphialign,minthealign,minpsialign);
				//
				i = -1;    
				itr = -1;
				syncro = 1;
				pktnum = 0;
				obt = 0;
				//
				while ( i < (nevents-1) ){    
				    //
				    // look for tracker data
				    //
				    Int_t ea = npia;
				    itr++;
				    i++;
				    syncro = 1;
				    if ( goodev[i] == 0 ) {
					goto skipthisev;					
				    };
				    ee = 0;
				    the = minthealign;
				    psi = minpsialign;
				    phi = minphialign;
				trkcalosync: 
				    //
				    // check if we have tracker data
				    //
				    if ( i >= nevents ){
					printf(" WARNING: no more calorimeter data.\n");
					goto closeandexit;
				    };
				    if ( itr > trnevents ){
					printf(" WARNING: no more tracker data.\n");
					goto closeandexit;
				    };
				    //
				    // retrieve tracker informations
				    //
				    tr->GetEntry(itr);
				    //
				    // check synchronization tracker and calorimeter informations:
				    //
				    otr->GetEntry(i);
				    //
				    pktnum = calo.pkt_num;
				    obt = calo.OBT;
				    if ( pktnum != trk.pkt_num || obt != trk.obt ){
					if ( pktnum > trk.pkt_num || obt > trk.obt ){
					    itr++;
					    syncro = 0;		
					    goto trkcalosync;
					};
					if ( pktnum < trk.pkt_num || obt < trk.obt ){
					    i++;
					    syncro = 0;
					    goto trkcalosync;
					};
				    };	
				    //
				    // here we have synchronized data
				    //
				    if ( !syncro ) {
					syncro = 1;
				    };       
				    //
				    if ( goodev[i] == 0 ) {
					goto skipthisev;					
				    };
				    ea = npia;
				    //
				    xcbar = (calo.cbar[ea][0] + minxalign) * cos(phi) * cos(the) + (calo.cbar[ea][1] + minyalign) * (cos(phi) * sin(the) * sin(psi) - sin(phi) * cos(psi)) + (piano[ea]  - zoff + minzalign) * (cos(phi) * sin(the) * cos(psi) + sin(phi) * sin (psi));
				    ycbar = (calo.cbar[ea][0] + minxalign) * sin(phi) * cos(the) + (calo.cbar[ea][1] + minyalign) * (sin(phi) * sin(the) * sin(phi) + cos(the) * cos(psi) ) + (piano[ea] - 5.1 - zoff + minzalign) * (sin(phi) * sin(the) * cos(psi) - cos(the) * sin(psi));
				    //
				    xtbar = mxtb[nza][i];
				    ytbar = -mytb[nza][i];
				    //
				    //
				    chi2nev += 1.;
				    //
				    chi2 += ((xcbar-xtbar)*(xcbar-xtbar)+(ycbar-ytbar)*(ycbar-ytbar))/(sigmacalo*sigmacalo);
				    //
				skipthisev:
				    continue;
				};
			    closeandexit:
				//
				minnxa = nxa;
				minnya = nya;
				minnza = nza;
				minphi = nphi;
				minpsi = npsi;
				minthe = nthe;
				tree->Fill();
				chi2nev = 0.;
				chi2 = 0.;
			    };
			};
		    };
		};
	    };
	};
	//
	// run over different alignment parameters
	//
	printf("\n\n");
	dxalig -= 7.;
	dyalig -= 7.;
	dzalig -= 7.;
	dphi -= 0.02;
	dpsi -= 0.02;
	dthe -= 0.02;
    };    
    hfile->Write();
    hfile->Close();
};


void findminimum(TString filename){
    const char* startingdir = gSystem->WorkingDirectory();
    TString path;
    stringcopy(path,startingdir);
    Int_t minnxa;
    Int_t minnza;
    Int_t minnya;
    Int_t minphi;
    Int_t minpsi;
    Int_t minthe;
    Float_t minxalign;
    Float_t minyalign;
    Float_t minzalign;
    Float_t minpsialign;
    Float_t minthealign;
    Float_t minphialign;
    Int_t nxa = 0;
    Int_t nza = 0;
    Int_t nya = 0;
    Int_t phi = 0;
    Int_t psi = 0;
    Int_t the = 0;
    Float_t chi2;
    Float_t chi2nev;
    TFile *f = new TFile(filename);
    TTree *tr = (TTree*) f->Get("CaloTrkAlign");
    tr->SetBranchAddress("chi2",&chi2);
    tr->SetBranchAddress("chi2nev",&chi2nev);
    tr->SetBranchAddress("minnxa",&minnxa);
    tr->SetBranchAddress("minnya",&minnya);
    tr->SetBranchAddress("minnza",&minnza);
    tr->SetBranchAddress("minphi",&minphi);
    tr->SetBranchAddress("minpsi",&minpsi);
    tr->SetBranchAddress("minthe",&minthe);
    tr->SetBranchAddress("minxalign",&minxalign);
    tr->SetBranchAddress("minyalign",&minyalign);
    tr->SetBranchAddress("minzalign",&minzalign);
    tr->SetBranchAddress("minpsialign",&minpsialign);
    tr->SetBranchAddress("minphialign",&minphialign);
    tr->SetBranchAddress("minthealign",&minthealign);
    //
    Long64_t nevents = tr->GetEntries();
    //    Float_t min = 1000000000.;
    Float_t min2 = 1000000000.;
    for (Int_t i = nevents-1; i>=0; i--){
	tr->GetEntry(i);
	//	if ( chi2/chi2nev <= min ){
	    //	    min = chi2/chi2nev;
	//    printf(" New minimum found:\n chi2 = %f minnxa = %i minnya = %i minnza = %i minthe = %i minpsi = %i minphi = %i \n\n",chi2/chi2nev,minnxa,minnya,minnza,minthe,minpsi,minphi); 	    
	//	};
	//	if ( minnxa==4 && minnya==6 && minnza==6 && chi2/chi2nev <= min2 ){
	if ( minnxa<6 && minnxa>2 && minnya<8 && minnya>4 && chi2/chi2nev <= min2 ){
	    min2 = chi2/chi2nev;
	    //	    printf(" 2- New minimum found:\n chi2 = %f minnxa = %i minnya = %i minnza = %i minthe = %i minpsi = %i minphi = %i \n\n",chi2/chi2nev,minnxa,minnya,minnza,minthe,minpsi,minphi); 	    
	    nxa = minnxa;
	    nya = minnya;
	    nza = minnza;
	    the = minthe;
	    phi = minphi;
	    psi = minpsi;
	};
    };
    printf("\n The best minimum found is for chi2 = %f with:\n  minnxa = %i minnya = %i minnza = %i minthe = %i minpsi = %i minphi = %i \n\n",min2,nxa,nya,nza,the,psi,phi);  
    Int_t nps = 0;
    Int_t nph = 0;
    Int_t nth = 0;
    Int_t nx = 0;
    Int_t ny = 0;
    Int_t nz = 0;
    Float_t vnpsi[11];
    Float_t vnphi[11];
    Float_t vnthe[11];
    Float_t vnx[11];
    Float_t vny[11];
    Float_t vnz[11];
    Float_t vrpsi[11];
    Float_t vrphi[11];
    Float_t vrthe[11];
    Float_t vrx[11];
    Float_t vry[11];
    Float_t vrz[11];
    for (Int_t i = nevents-1; i>=0; i--){
	tr->GetEntry(i);    
	if ( minnxa==nxa && minnya==nya && minnza==nza && minthe==the && minphi==phi ){
	    vnpsi[nps] = chi2/chi2nev;
	    vrpsi[nps] = minpsialign;
	    nps++;
	};
	if ( minnxa==nxa && minnya==nya && minnza==nza && minthe==the && minpsi==psi ){
	    vnphi[nph] = chi2/chi2nev;
	    vrphi[nph] = minphialign;
	    nph++;
	};
	if ( minnxa==nxa && minnya==nya && minnza==nza && minphi==phi && minpsi==psi ){
	    vnthe[nth] = chi2/chi2nev;
	    vrthe[nth] = minthealign;
	    nth++;
	};
	if ( minnxa==nxa && minnya==nya && minthe==the && minphi==phi && minpsi==psi ){
	    vnz[nz] = chi2/chi2nev;
	    vrz[nz] = minzalign;
	    nz++;
	};
	if ( minnxa==nxa && minnza==nza && minthe==the && minphi==phi && minpsi==psi ){
	    vny[ny] = chi2/chi2nev;
	    vry[ny] = -minyalign;
	    ny++;
	};
	if ( minnya==nya && minnza==nza && minthe==the && minphi==phi && minpsi==psi ){
	    vnx[nx] = chi2/chi2nev;
	    vrx[nx] = -minxalign;
	    nx++;
	};
    };
    //
    //
    //
    const string myfil = (const char *)filename;
    Int_t myposiz = myfil.find("dw_");    
    if ( myposiz == -1 ) {
      myposiz = myfil.find("DW_");
    };
    TString fileid;
    stringcopy(fileid,filename,myposiz+18,myposiz+28);
    //
    //
    stringstream namefile;
    const char *pattho = path;
    //
    TCanvas *figura1 = new TCanvas("figura1","figura1", 1100, 900);
    TGraph *gr;
    TF1 *fit;
    Float_t par2;
    Float_t par3;
    stringstream titolo;
    //
    figura1->Clear();
    figura1->cd();
    gPad->SetTickx();
    gPad->SetTicky();
    gr = new TGraph(11,vrx,vnx);
    gr->SetLineWidth(1);
    gr->SetMarkerStyle(20);
    gStyle->SetOptFit(0);
    gr->GetXaxis()->SetTitle("x offset [mm]");
    gr->GetYaxis()->SetTitle("#chi^{2}");
    gr->Fit("pol2","q");
    fit = gr->GetFunction("pol2"); 
    par2 = (float)fit->GetParameter(1);
    par3 = (float)fit->GetParameter(2);
    titolo.str("");
    titolo << " Minimum at x = " << -par2/(2.*par3);
    titolo << " mm";
    gr->SetTitle(titolo.str().c_str());
    gr->Draw("AP");    
    namefile.str("");
    namefile << pattho << "/x_chi2_";
    namefile << fileid.Data()  << ".ps";
    printf(" Save figure %s\n",namefile.str().c_str());
    figura1->SaveAs(namefile.str().c_str());
    //
    figura1->Clear();
    figura1->cd();
    gPad->SetTickx();
    gPad->SetTicky();
    gr = new TGraph(11,vry,vny);
    gr->SetLineWidth(1);
    gr->SetMarkerStyle(20);
    gStyle->SetOptFit(0);
    gr->GetXaxis()->SetTitle("y offset [mm]");
    gr->GetYaxis()->SetTitle("#chi^{2}");
    gr->Fit("pol2","q");
    fit = gr->GetFunction("pol2"); 
    par2 = (float)fit->GetParameter(1);
    par3 = (float)fit->GetParameter(2);
    titolo.str("");
    titolo << " Minimum at y = " << -par2/(2.*par3);
    titolo << " mm";
    gr->SetTitle(titolo.str().c_str());
    gr->Draw("AP");    
    namefile.str("");
    namefile << pattho << "/y_chi2_";
    namefile << fileid.Data()  << ".ps";
    printf(" Save figure %s\n",namefile.str().c_str());
    figura1->SaveAs(namefile.str().c_str());
    //
    figura1->Clear();
    figura1->cd();
    gPad->SetTickx();
    gPad->SetTicky();
    gr = new TGraph(11,vrz,vnz);
    gr->SetLineWidth(1);
    gr->SetMarkerStyle(20);
    gStyle->SetOptFit(0);
    gr->GetXaxis()->SetTitle("z offset [mm]");
    gr->GetYaxis()->SetTitle("#chi^{2}");
    gr->Fit("pol2","q");
    fit = gr->GetFunction("pol2"); 
    par2 = (float)fit->GetParameter(1);
    par3 = (float)fit->GetParameter(2);
    titolo.str("");
    titolo << " Minimum at z = " << -par2/(2.*par3);
    titolo << " mm";
    gr->SetTitle(titolo.str().c_str());
    gr->Draw("AP");    
    namefile.str("");
    namefile << pattho << "/z_chi2_";
    namefile << fileid.Data()  << ".ps";
    printf(" Save figure %s\n",namefile.str().c_str());
    figura1->SaveAs(namefile.str().c_str());
    //
    figura1->Clear();
    figura1->cd();
    gPad->SetTickx();
    gPad->SetTicky();
    gr = new TGraph(11,vrphi,vnphi);
    gr->SetLineWidth(1);
    gr->SetMarkerStyle(20);
    gStyle->SetOptFit(0);
    gr->GetXaxis()->SetTitle("#phi_{z} offset [rad]");
    gr->GetYaxis()->SetTitle("#chi^{2}");
    gr->Fit("pol2","q");
    fit = gr->GetFunction("pol2"); 
    par2 = (float)fit->GetParameter(1);
    par3 = (float)fit->GetParameter(2);
    titolo.str("");
    titolo << " Minimum at #phi_{z} = " << -par2/(2.*par3);
    titolo << " rad";
    gr->SetTitle(titolo.str().c_str());
    gr->Draw("AP");    
    namefile.str("");
    namefile << pattho << "/phi_chi2_";
    namefile << fileid.Data()  << ".ps";
    printf(" Save figure %s\n",namefile.str().c_str());
    figura1->SaveAs(namefile.str().c_str());
    //
    figura1->Clear();
    figura1->cd();
    gPad->SetTickx();
    gPad->SetTicky();
    gr = new TGraph(11,vrthe,vnthe);
    gr->SetLineWidth(1);
    gr->SetMarkerStyle(20);
    gStyle->SetOptFit(0);
    gr->GetXaxis()->SetTitle("#theta_{y} offset [rad]");
    gr->GetYaxis()->SetTitle("#chi^{2}");
    gr->Fit("pol2","q");
    fit = gr->GetFunction("pol2"); 
    par2 = (float)fit->GetParameter(1);
    par3 = (float)fit->GetParameter(2);
    titolo.str("");
    titolo << " Minimum at #theta_{y} = " << -par2/(2.*par3);
    titolo << " rad";
    gr->SetTitle(titolo.str().c_str());
    gr->Draw("AP");    
    namefile.str("");
    namefile << pattho << "/the_chi2_";
    namefile << fileid.Data()  << ".ps";
    printf(" Save figure %s\n",namefile.str().c_str());
    figura1->SaveAs(namefile.str().c_str());
    //
    figura1->Clear();
    figura1->cd();
    gPad->SetTickx();
    gPad->SetTicky();
    gr = new TGraph(11,vrpsi,vnpsi);
    gr->SetLineWidth(1);
    gr->SetMarkerStyle(20);
    gStyle->SetOptFit(0);
    gr->GetXaxis()->SetTitle("#psi_{x} offset [rad]");
    gr->GetYaxis()->SetTitle("#chi^{2}");
    gr->Fit("pol2","q");
    fit = gr->GetFunction("pol2"); 
    par2 = (float)fit->GetParameter(1);
    par3 = (float)fit->GetParameter(2);
    titolo.str("");
    titolo << " Minimum at #psi_{x} = " << -par2/(2.*par3);
    titolo << " rad";
    gr->SetTitle(titolo.str().c_str());
    gr->Draw("AP");
    namefile.str("");
    namefile << pattho << "/psi_chi2_";
    namefile << fileid.Data()  << ".ps";
    printf(" Save figure %s\n",namefile.str().c_str());
    figura1->SaveAs(namefile.str().c_str());
    //
}

void findminimumnorotation(TString filename){
    Int_t minnxa;
    Int_t minnza;
    Int_t minnya;
    Int_t minphi;
    Int_t minpsi;
    Int_t minthe;
    Float_t chi2;
    Float_t chi2nev;
    TFile *f = new TFile(filename);
    TTree *tr = (TTree*) f->Get("CaloTrkAlign");
    tr->SetBranchAddress("chi2",&chi2);
    tr->SetBranchAddress("chi2nev",&chi2nev);
    tr->SetBranchAddress("minnxa",&minnxa);
    tr->SetBranchAddress("minnya",&minnya);
    tr->SetBranchAddress("minnza",&minnza);
    tr->SetBranchAddress("minphi",&minphi);
    tr->SetBranchAddress("minpsi",&minpsi);
    tr->SetBranchAddress("minthe",&minthe);
    Long64_t nevents = tr->GetEntries();
    Float_t min = 1000000000.;
    for (Int_t i = nevents-1; i>=0; i--){
	tr->GetEntry(i);
	if ( chi2/chi2nev <= min && minthe==5 && minpsi==5 && minphi == 5 ){
	    min = chi2/chi2nev;
	    printf(" New minimum found:\n chi2 = %f minnxa = %i minnya = %i minnza = %i minthe = %i minpsi = %i minphi = %i \n\n",chi2/chi2nev,minnxa,minnya,minnza,minthe,minpsi,minphi); 	    
	};
    };
    

}