//Plotting macro for data. 
//Needs to be compiled!
//Author: Tiger

#include "TCanvas.h"
#include "TFile.h"
#include "TH2F.h"
#include "TH3F.h"
#include "TClonesArray.h"
#include "TGeoManager.h"
#include "TTree.h"
#include "TObjArray.h"
#include "TPolyMarker.h"
#include "TPolyMarker3D.h"
#include "TPolyLine.h"
#include "TVector3.h"
#include "TBox.h"
#include "TSystem.h"
#include "TROOT.h"
#include "TMinuit.h"
#include "tpc/par/TpcGem.h"
#include "tpc/par/TpcPadPlane.h"
#include "tpc/par/TpcPadShapePool.h"
#include "tpc/base/TpcDigi.h"
#include "tpc/par/TpcDigiMapper.h"
#include "tpc/par/TpcGas.h"
#include "TView3D.h"
#include "TApplication.h"
#include "TMath.h"
#include "TStyle.h"
#include "TString.h"
#include "TColor.h"
#include "THStack.h"
#include "TVector3.h"
#include "tpc/par/TpcAlignmentManager.h"
//#include "tpc/TestBench/TBStripCluster.h"
#include "tpc/base/TpcCluster.h"
#include "tpc/TestBench/TBGEMCluster.h"
#include "tpc/base/TpcResidual.h"
#include "tpc/TestBench/TBSICluster.h"
#include "tpc/tools/residuals/TpcRefResidualCollection.h"
#include "tpc/base/TpcResidualCollection.h"
#include "tpc/tools/residuals/TpcRefGFTrkResCalc.h"
#include "GFAbsTrackRep.h"
#include "GFAbsRecoHit.h"
#include "GFConstField.h"
#include "GFDetPlane.h"
#include "GFException.h"
#include "GFFieldManager.h"
#include "GFKalman.h"
#include "GFTools.h"
#include "GFTrack.h"
#include "GFTrackCand.h"
//#include "GenfitTools/trackrep/RKTrackRep/RKTrackRep.h"


#include <vector>
#include <cmath>
#include <iostream>
#include <algorithm>
#include <map>
#include <string>
#include <exception>

#define DEBUG 1
const int NDET = 8;
const int NDETXY = 4;
const int TPCDET = 8;
const int NSIZE = 11;
//class TBTSICluster;
//GLOBALS FOR MINIMAZATION
std::vector<TpcResidual*> vres;
TpcAlignmentManager * alMan;
Double_t diffufit(Double_t * x , Double_t*par);

//
void plotTCres(TString filename,  unsigned int EvMax=0,std::string alfile = "macro/tpc/TestBench/COMPASS_alignment.tpcAligned.txt",std::string trackarray="TBTrackPostFit")
{
  TGeoManager* geom = new TGeoManager("Geometry", "Geane geometry");
  TGeoManager::Import("tpc/parfiles/FOPIGeo.root");
  int lol = 0;
  std::cout << "lol " << lol << std::endl; lol++;
  GFFieldManager::getInstance()->init(new GFConstField(0.,0.,0.));
  std::cout << "lol " << lol << std::endl; lol++;
  while(gROOT->GetListOfCanvases()->GetSize())
    delete gROOT->GetListOfCanvases()->At(0);
  gStyle->SetOptStat(10);
  gStyle->SetPalette(1);
  int a = 0;
  TFile* recofile = new TFile(filename);
  if(recofile->IsZombie()) {
    std::cerr<<"Reco file not existing! Aborting."<<std::endl;
    return;  }
  std::cout << "lol " << lol << std::endl; lol++;
  TTree* recotree = (TTree*)recofile->Get("cbmsim");
  if(!recotree) {
    std::cerr<<"Reco Tree not existing! Aborting."<<std::endl;
    return;  }

 std::map<TString, TClonesArray*> Tbranches;
 
 //I CAN HAZ ARRAYS :3
 // TClonesArray* trTBArr=new TClonesArray("GFTrack");
 // recotree->SetBranchAddress("TBTrackPostFit", &trTBArr);
 // Tbranches["TBTrackPostFit"] = trTBArr;
 // if(trTBArr==0) {std::cout<< "Init TB GFTrack-Array not found!"<< std::endl;return;}

 TClonesArray* trTBArr=new TClonesArray("GFTrack");
 recotree->SetBranchAddress(trackarray.c_str(), &trTBArr);
 Tbranches[trackarray.c_str()] = trTBArr;
 if(trTBArr==0) {std::cout<< "Init TB GFTrack-Array not found!"<< std::endl;return;}
 
 TClonesArray* tpcArr=new TClonesArray("TpcCluster");
 recotree->SetBranchAddress("TpcCluster", &tpcArr);
 Tbranches["TpcCluster"] = tpcArr;
 if(tpcArr==0) {std::cout<< "Init TpcCluster Array not found!"<< std::endl;return;}

 TClonesArray* tpcResArr=new TClonesArray("TpcRefResidualCollection");
 recotree->SetBranchAddress("TBTPCres", &tpcResArr);
 if(tpcResArr==0) {std::cout<< "Init TpcResidualCollection Array not found!"<< std::endl;return;}
 std::cout << "Loading alMan " << std::endl; 
 alMan = TpcAlignmentManager::getInstance(alfile.c_str());
 if(alMan==0) {std::cout<< "Init TpcAlignmentManager:"<< alfile.c_str()<<" not found!"<< std::endl;return;}
 std::cout << "Turing test passed." << std::endl;
 int nbinmap = 300;
   
  //TH1...
 TH1F * hresx = new TH1F("hresx","hresx",400,-0.2,0.2);
 TH1F * hresy = new TH1F("hresy","hresy",400,-0.2,0.2);
 TH1F * hresz = new TH1F("hresz","hresz",400,-0.3,0.3);
 TH1F * hresdt = new TH1F("hresdt","hresdt",800,-3500,0);
 TH1F * hrescldt = new TH1F("hrescldt","hrescldt",800,-3500,0);
  char blamiaou[500];std::vector <TH1F*> vhresx, vhresy, vhresz;
  for(int o = 1;o<5;o++)
    {
      sprintf(blamiaou,"hresx_Size%d",o);
      TH1F * hresxx = new TH1F(blamiaou,blamiaou,400,-0.2,0.2);
      sprintf(blamiaou,"hresy_Size%d",o);
      TH1F * hresyy = new TH1F(blamiaou,blamiaou,400,-0.2,0.2);
      sprintf(blamiaou,"hresz_Size%d",o);
      TH1F * hreszz = new TH1F(blamiaou,blamiaou,400,-0.3,0.3);
      switch(o)
	{
	case 1 : {hresxx->SetLineColor(kRed); hresyy->SetLineColor(kRed); hreszz->SetLineColor(kRed);} break;
	case 2 : {hresxx->SetLineColor(kGreen+4); hresyy->SetLineColor(kGreen+4); hreszz->SetLineColor(kGreen+4);} break;
	case 3 : {hresxx->SetLineColor(kCyan+3); hresyy->SetLineColor(kCyan+3); hreszz->SetLineColor(kCyan+3);} break;
	case 4 : {hresxx->SetLineColor(kBlue); hresyy->SetLineColor(kBlue); hreszz->SetLineColor(kBlue);} break;	}
      vhresx.push_back(hresxx);
      vhresy.push_back(hresyy);
      vhresz.push_back(hreszz);
      std::cout << blamiaou << std::endl;    }

  TH2F * hresxy = new TH2F("hresxy","hresxy", 800,-0.15,0.15,800,-0.15,0.15);
  TH2F * hresyz = new TH2F("hresyz","hresyz", 800,-0.15,0.15,800,-0.15,0.15);
  
  TH2F * hresxvsx = new TH2F("hresxvsx","hresxvsx", 300, -4.7,4.7,300,-0.1,0.1);
  TH2F * hresxvsz = new TH2F("hresxvsz","hresxvsz", 300, -4.7,4.7,300,-0.1,0.1);
  TH2F * hresyvsy = new TH2F("hresyvsy","hresyvsy", 300, 0.7,2.7,300,-0.1,0.1);
  TH2F * hreszvsz = new TH2F("hreszvsz","hreszvsz", 25, 2.5,7.,300,-0.1,0.1);
  TH2F * hresyvsz = new TH2F("hresyvsz","hresyvsz", 25, 2.5,7.,300,-0.1,0.1);
    
  TH2F * hampratio = new TH2F("hampratio","hampratio",30,0,30,500,-0.05,1.05);


  std::vector <TH2F*> vhpadresp;
  for(int o = 1;o<NSIZE+2;o++)
    {
      sprintf(blamiaou,"hpadresp_Size%d",o);
      TH2F * hpadresp = new TH2F(blamiaou,blamiaou,150,-0.15,0.15,150,-0.15,0.15);
      vhpadresp.push_back(hpadresp);
    }
  TH2F * hitmap = new TH2F("hitmap","hitmap",nbinmap,-4.7,4.7,nbinmap,-4.7,4.7);
  TH2F * hitmapxz = new TH2F("hitmapxz","hitmapxz",nbinmap,-4.7,4.7,nbinmap,0,10);
  TH2F * hittrmap = new TH2F("hittrmap","hittrmap",nbinmap,-4.7,4.7,nbinmap,-4.7,4.7);
  TH2F * hittrmapxz = new TH2F("hittrmapxz","hittrmapxz",nbinmap,-4.7,4.7,nbinmap,0,10);
  TH2F * hitmapmag = new TH2F("hitmapmag","hitmapmag",nbinmap,-4.7,4.7,nbinmap,-4.7,4.7);
  TH2F * hitmapxzmag = new TH2F("hitmapxzmag","hitmapxzmag",nbinmap,-4.7,4.7,nbinmap,0,10);
  TH3F * hitmapxyz = new TH3F("hitmapxyz","hitmapxyz",100,-4.7,4.7,100,-4.7,-4.7,100,0,10);
  //Residuals calculator
  
  TpcRefGFTrkResCalc* resgft = new TpcRefGFTrkResCalc();
  int test =   resgft->addBranchName("TpcCluster");
  if(test) {std::cout << "resgft->addBranchName(\"TpcCluster\"); is inconsistent" << std::endl;}
  //  test = resgft->addBranchName("TBTrackPostFit");
  test = resgft->addBranchName(trackarray.c_str());
  //  if(test) {std::cout << "resgft->addBranchName(\"TBTrackPostFit\"); is inconsistent" << std::endl;}
  if(test) {std::cout << "resgft->addBranchName(\"TrackPostFit\"); is inconsistent" << std::endl;}
  TVector3 uu(1.,0.,0.);  TVector3 vv(0.,1.,0.);  TVector3 oo(0.,0.,105.);//plane in front of the TC
  resgft->setProjectionPlane(oo,uu,vv);
  resgft->setTrackRepId(0);
  std::cout << "lol " << lol << std::endl; lol++;
  std::map<TString, TClonesArray*>* branches = resgft->getBranchMap();
  std::map<TString, TClonesArray*>::iterator it;
  std::cout << "lo---lul " << lol << std::endl; lol++;
  for(it=branches->begin(); it!=branches->end(); it++) {
    TClonesArray* br = Tbranches[it->first];
    if(br==NULL) {
      std::string s("Branch ");
      s.append(it->first.Data());
      s.append(" requested by residual calculator, but not found in tree!");
      Fatal("TpcRefTrackResidualTask::Init()",s.c_str());
    }
    it->second = br;
  } 

  std::cout << "digi mapper innit : " << std::endl;
  const TpcGem* _gem=new TpcGem(5000, 0.02);          // Gain, Spread
  const TpcPadShapePool* _padShapes = new TpcPadShapePool("tpc/parfiles/TBhexa_pads.dat",  *_gem,  0.5, 0.02,  0.01); 
  TpcPadPlane* _padPlane= new TpcPadPlane("tpc/parfiles/newPadTC.dat",_padShapes);
  const TpcGas* gas=new TpcGas("tpc/parfiles/ARGON-70_CO2-30_B0.1_PRES1013.asc",400);
  TpcDigiMapper::getInstance(false)->init(_padPlane,_gem,gas, _padShapes,-40.,-60000.,40.);
  
  std::cout << "Self awareness ... " << std::endl;
  resgft->init();
  unsigned int MAX =recotree->GetEntries();
  std::cout << "Number of events: " << MAX  <<std::endl;
  if(!EvMax)EvMax=MAX;
  vres.clear();
  for(unsigned int k=0; k<MAX && k<EvMax; k++) 
    {
      recotree->GetEvent(k);
      int ntpccl = tpcArr->GetEntries();
      resgft->calc();
      std::vector<TpcRefResidualCollection*>* blka =  resgft->getResiduals() ;
      int nRes = blka->size();
      for(int ires=0; ires<nRes; ires++) 
	{
	  TpcRefResidualCollection * resCol = (*blka)[ires];
	  for (int ii = 0 ; ii < resCol->getSize(); ii++)
	    {
	      TpcResidual * miaou = new TpcResidual(*resCol->getResidual(ii));
	      TpcCluster* cl = (TpcCluster*)tpcArr->At(miaou->getHitIndex());
	      
	      //	      if(cl->amp()<1000) continue;//cut on bg
	      //vres.push_back(miaou);
	      TVector3 miaouv = alMan->masterToLocalVect("TC01____",miaou->getResXYZ());
	      TVector3 miaoupos = alMan->masterToLocal("TC01____",miaou->getHitPos());
	      TVector3 miaourefpos = alMan->masterToLocal("TC01____",miaou->getRefPos());
	      double corrz = miaoupos.Z()*-1.55e-2+7.32000000000000012e-02;
	      miaoupos.SetZ(miaoupos.Z()-corrz);

	      hitmap->Fill(miaoupos.X(),miaoupos.Y());
	      hitmapxz->Fill(miaoupos.X(),miaoupos.Z());
	      hittrmap->Fill(miaourefpos.X(),miaourefpos.Y());
	      hittrmapxz->Fill(miaourefpos.X(),miaourefpos.Z());
	      hitmapxyz->Fill(miaourefpos.X(),miaourefpos.Y(),miaourefpos.Z());
	      hitmapmag->Fill(miaoupos.X(),miaoupos.Y(),miaouv.Mag());
	      hitmapxzmag->Fill(miaoupos.X(),miaoupos.Z(),miaouv.Mag());
	      //small pads

	      if(1)//miaoupos.Y()>0.6)
		{//Y U NO CUT
		//if(miaoupos.Y()<0.6){
		
		int clsi = cl->size(); 
		hresx->Fill(miaouv.X());
		hresy->Fill(miaouv.Y());

		if(miaoupos.Z()>3.5&&miaoupos.Z()<4)
		  hresz->Fill(miaoupos.Z()-miaourefpos.Z());//miaouv.Z());

		hresxy->Fill(miaouv.X(),miaouv.Y());
		hresyz->Fill(miaouv.Y(),miaoupos.Z()-miaourefpos.Z());//miaouv.Z());

		hresxvsx->Fill(miaoupos.X(),miaouv.X());
		hresxvsz->Fill(miaoupos.X(),miaoupos.Z()-miaourefpos.Z());//miaouv.Z());
		hresyvsy->Fill(miaoupos.Y(),miaouv.Y());
		//		hreszvsz->Fill(miaoupos.Z(),miaouv.Z());
		hreszvsz->Fill(miaoupos.Z(),miaoupos.Z()-miaourefpos.Z());
		hresyvsz->Fill(miaoupos.Z(),miaoupos.Y()-miaourefpos.Y());

		if(clsi>0)
		  if(clsi<4)
		    {
		      vhresx.at(clsi-1)->Fill(miaouv.X());
		      vhresy.at(clsi-1)->Fill(miaouv.Y());
		      vhresz.at(clsi-1)->Fill(miaoupos.Z()-miaourefpos.Z());//miaouv.Z());
		    }
		  else 
		    {
		      vhresx.at(3)->Fill(miaouv.X());
		      vhresy.at(3)->Fill(miaouv.Y());
		      vhresz.at(3)->Fill(miaoupos.Z()-miaourefpos.Z());//miaouv.Z());
		    }

		//corresponding cluster
		TVector3 digpos; double maxampdigi = 0;
		double tcm = 0.00055192;//0.00055 cm/ns
		for (int mew = 0; mew < clsi ; mew++)
		  {//digi loop
		  hampratio->Fill(cl->size(),cl->getDigi(mew)->amp()/cl->amp());
		  if(cl->posuvw().Z()>3.2 &&cl->posuvw().Z()<3.8)//slice in z to avoid calib issue
		    if(cl->getDigi(mew)->tlength()>5)
		      {
			hresdt->Fill(miaourefpos.Z()/tcm-cl->getDigi(mew)->t()*60.);
			hrescldt->Fill(miaoupos.Z()/tcm-cl->getDigi(mew)->t()*60.);
		      }
		  if(cl->getDigi(mew)->amp()>maxampdigi)
		    {
		      maxampdigi = cl->getDigi(mew)->amp();
		      TpcDigiMapper::getInstance()->map(cl->getDigi(mew),digpos);
		    }}
		TVector3 padrespv = digpos-cl->posuvw();
		if(clsi>0 && clsi<NSIZE)
		  vhpadresp.at(clsi)->Fill(padrespv.X(),padrespv.Y());
		if(clsi>0)
		  vhpadresp.at(NSIZE)->Fill(padrespv.X(),padrespv.Y());
	      }
	    }
	}
    }//end event loop 
  std::cout << "DRAWING " << std::endl;
  
  //Drawing() 
  
  gROOT->SetStyle("Plain");
  gStyle->SetPalette(1);
  const Int_t NCont = 250;
  const Int_t NRGBs = 5;
  Double_t stops[NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 };
  Double_t red[NRGBs]   = { 0.00, 0.00, 0.87, 1.00, 0.51 };
  Double_t green[NRGBs] = { 0.00, 0.81, 1.00, 0.20, 0.00 };
  Double_t blue[NRGBs]  = { 0.51, 1.00, 0.12, 0.00, 0.00 };
  TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);

  TCanvas * cx = new TCanvas("cx","cx",1200,600);
  cx->Divide(3);
  cx->cd(1);hresx->Draw();  for(int o = 0;o<4;o++) vhresx.at(o)->Draw("same");
  cx->cd(2);hresy->Draw();  for(int o = 0;o<4;o++) vhresy.at(o)->Draw("same");
  cx->cd(3);hresz->Draw();  for(int o = 0;o<4;o++) vhresz.at(o)->Draw("same");
  
  TCanvas * cxvsx = new TCanvas("cxvsx","cxvsx",1200,600);
  cxvsx->Divide(3);
  cxvsx->cd(1);hresxvsx->Draw("colz");//hresxvsx->Draw("colz");
  cxvsx->cd(2);hresyvsy->Draw("colz");
  cxvsx->cd(3);hreszvsz->Draw("colz");

  TCanvas * csl = new TCanvas("csl","Diffusion",1200,600);
  csl->Divide(1,2);
  csl->cd(1);
  TObjArray slices;
  hreszvsz->FitSlicesY(0,0,-1,0,"QNR",&slices);
  TH1D * hsl = (TH1D*)slices.At(2);//(TH1D*)gDirectory->GetObject("hreszvsz_2");
  hsl->DrawCopy();
  csl->cd(2);
  TObjArray slices2;
  hresyvsz->FitSlicesY(0,0,-1,0,"QNR",&slices2);
  TH1D * hsl2 = (TH1D*)slices2.At(2);//(TH1D*)gDirectory->GetObject("hreszvsz_2");
  TH1D * hsl2_p = (TH1D*)hsl2->Clone("hsl2_p");
  hsl2_p->Draw();
  TF1 * f2 = new TF1("sqrtFit",diffufit,1,10,4);
  f2->SetLineColor(kBlue);
  f2->SetParNames("reso","Dt","Neff","zoff");
  f2->SetParLimits(0,0,0.099);
  f2->SetParLimits(1,0,0.5);
  f2->SetParLimits(2,0.,500.);
  f2->SetParLimits(3,-10,10.);
  f2->SetNpx(5000); 
  TFitResultPtr r = hsl2_p->Fit(f2,"","",3,7);
	      
  TCanvas * cxy = new TCanvas("cxy","cxy",1200,600);
  cxy->Divide(2);
  cxy->cd(1);hresxy->Draw("colz");
  cxy->cd(2);hresyz->Draw("colz");
  
  TCanvas * cpadresp = new TCanvas("cpadresp","cpadresp",1200,800);
  cpadresp->Divide((NSIZE+1)/2,2); 
  for(int o = 1; o< NSIZE+1;o++)
    {
      cpadresp->cd(o);
      std::cout << "ok = "<< o << std::endl;
      vhpadresp.at(o)->Draw("colz");
      std::cout << "ok = "<< o << std::endl;
    }

  TCanvas * cmap = new TCanvas("cmap","cmap",1200,800);
  cmap->Divide(2);
  cmap->cd(1);hitmap->Draw("colz");
  cmap->cd(2);hitmapxz->Draw("colz");

  TCanvas * cxyz = new TCanvas("cxyz","cxyz",1200,800);
  hitmapxyz->GetXaxis()->SetRangeUser(-4.7,4.7);
  hitmapxyz->GetYaxis()->SetRangeUser(-4.7,4.7);
  hitmapxyz->GetZaxis()->SetRangeUser(0,10);
  hitmapxyz->Draw();

  TCanvas * ctrmap = new TCanvas("ctrmap","ctrmap",1200,800);
  ctrmap->Divide(2);
  ctrmap->cd(1);hittrmap->Draw("colz");
  ctrmap->cd(2);hittrmapxz->Draw("colz");

  TCanvas * cdig = new TCanvas("digitime","digitime",1200,800);
  cdig->Divide(2);
  cdig->cd(1); hresdt->Draw();
  cdig->cd(2); hrescldt->Draw();


  TCanvas * crat = new TCanvas("crat","crat",1200,800);
  hampratio->Draw("colz");
  TCanvas * cmapmag = new TCanvas("cmapmag","cmapmag",1200,800);
  cmapmag->Divide(2);
  cmapmag->cd(1);hitmapmag->Divide(hitmap);hitmapmag->Draw("colz");
  cmapmag->cd(2);hitmapxzmag->Divide(hitmapxz);hitmapxzmag->Draw("colz");

  //residuals + FIT
  TCanvas * cxX = new TCanvas("cxX","cxX",1200,600);
  cxX->Divide(3);
  cxX->cd(1);hresx->Draw();
  TF1* fgpresx2 = new TF1("fgpresx2", "gaus(0)");
  fgpresx2->SetParName(0, "Ag");
  double max2 = hresx->GetMaximum();
  fgpresx2->SetParameter(0, max2);
  fgpresx2->SetParName(1, "mean1");
  fgpresx2->SetParameter(1, 0.);
  fgpresx2->SetParName(2, "sigma1");
  fgpresx2->SetParameter(2, 0.03);
  hresx->Fit("fgpresx2","","",-0.03,0.03);
  hresx->GetXaxis()->SetRangeUser(-0.06,0.06);
  
  TF1* fggpresx2 = new TF1("fggpresx2", "gaus(0) +  gaus(3)");
  fggpresx2->SetParName(0, "Ag");
  fggpresx2->SetParameter(0, fgpresx2->GetParameter(0));
  fggpresx2->SetParName(1, "mean1");
  fggpresx2->SetParameter(1, fgpresx2->GetParameter(1) );
  fggpresx2->SetParName(2, "sigma1");
  fggpresx2->SetParameter(2, fabs(fgpresx2->GetParameter(2)));
  fggpresx2->SetParName(3, "Ag2");
  fggpresx2->SetParameter(3, max2);
  fggpresx2->SetParName(4, "mean2");
  fggpresx2->SetParameter(4, 0.);
  fggpresx2->SetParName(5, "sigma2");
  fggpresx2->SetParameter(5, 0.1);
  fggpresx2->SetLineWidth(2); fggpresx2->SetNpx(1000);
  hresx->Fit("fggpresx2","","",-0.1,0.1);
  hresx->Draw();
  //draw gaussians
  TF1* fg1resx2 = new TF1("fg1resx2", "gaus", -0.2,0.2);
  for(int i = 0;i<3;++i) fg1resx2->SetParameter(i,fggpresx2->GetParameter(i));
  fg1resx2->SetLineColor(kBlue);fg1resx2->SetLineWidth(2); fg1resx2->SetNpx(1000);
  fg1resx2->Draw("same");
  
  TF1* fg2resx2 = new TF1("fg2resx2","gaus", -0.2,0.2);
  for(int i = 0;i<3;++i) fg2resx2->SetParameter(i,fggpresx2->GetParameter(i+3));
  fg2resx2->SetLineColor(kRed);fg2resx2->SetLineWidth(2); fg2resx2->SetNpx(1000);
  fg2resx2->Draw("same");
 
  double dresx2 = fabs(fggpresx2->GetParameter(2));double ampresx2 = fggpresx2->GetParameter(0);
  double dresx3 = fabs(fggpresx2->GetParameter(5));double ampresx3 = fggpresx2->GetParameter(3);
  std::cout << dresx2 << " x " << ampresx2 << "\t" << dresx3 << " x " << ampresx3 << endl;
  double RESresx = (dresx2*ampresx2+dresx3*ampresx3) / (ampresx2+ampresx3);	
  std::cout << "== " << (dresx2*ampresx2+dresx3*ampresx3) / (ampresx2+ampresx3) << std::endl;
  char titleresx2[560];sprintf(titleresx2,"residualU, #sigma = %f #mum",RESresx*10000.);
  hresx->SetTitle(titleresx2);
  hresx->GetXaxis()->SetTitle("Res. (cm)");
  //////////////////////////////////////////////////////////////////////////////////////////////////////////////
  cxX->cd(2);hresy->Draw();
  TF1* fgpresy2 = new TF1("fgpresy2", "gaus(0)");
  fgpresy2->SetParName(0, "Ag");
  double max23 = hresy->GetMaximum();
  fgpresy2->SetParameter(0, max23);
  fgpresy2->SetParName(1, "mean1");
  fgpresy2->SetParameter(1, 0.);
  fgpresy2->SetParName(2, "sigma1");
  fgpresy2->SetParameter(2, 0.02);
  hresy->Fit("fgpresy2","","",-0.03,0.03);
  hresy->GetXaxis()->SetRangeUser(-0.04,0.04);
  
  TF1* fggpresy2 = new TF1("fggpresy2", "gaus(0) +  gaus(3)");
  fggpresy2->SetParName(0, "Ag");
  fggpresy2->SetParameter(0, fgpresy2->GetParameter(0));
  fggpresy2->SetParName(1, "mean1");
  fggpresy2->SetParameter(1, fgpresy2->GetParameter(1) );
  fggpresy2->SetParName(2, "sigma1");
  fggpresy2->SetParameter(2, fabs(fgpresy2->GetParameter(2)));
  fggpresy2->SetParName(3, "Ag2");
  fggpresy2->SetParameter(3, max23);
  fggpresy2->SetParName(4, "mean2");
  fggpresy2->SetParameter(4, 0.);
  fggpresy2->SetParName(5, "sigma2");
  fggpresy2->SetParameter(5, 0.1);
  fggpresy2->SetLineWidth(2); fggpresy2->SetNpx(1000);
  hresy->GetXaxis()->SetRangeUser(-0.2,0.2);
  hresy->Fit("fggpresy2","","",-0.2,0.2);
  hresy->Draw();
  //draw gaussians
  TF1* fg1resy2 = new TF1("fg1resy2", "gaus", -0.2,0.2);
  for(int i = 0;i<3;++i) fg1resy2->SetParameter(i,fggpresy2->GetParameter(i));
  fg1resy2->SetLineColor(kBlue);fg1resy2->SetLineWidth(2); fg1resy2->SetNpx(1000);
  fg1resy2->Draw("same");
  
  TF1* fg2resy2 = new TF1("fg2resy2","gaus", -0.2,0.2);
  for(int i = 0;i<3;++i) fg2resy2->SetParameter(i,fggpresy2->GetParameter(i+3));
  fg2resy2->SetLineColor(kRed);fg2resy2->SetLineWidth(2); fg2resy2->SetNpx(1000);
  fg2resy2->Draw("same");
 
  double dresy2 = fabs(fggpresy2->GetParameter(2));double ampresy2 = fggpresy2->GetParameter(0);
  double dresy3 = fabs(fggpresy2->GetParameter(5));double ampresy3 = fggpresy2->GetParameter(3);
  std::cout << dresy2 << " x " << ampresy2 << "\t" << dresy3 << " x " << ampresy3 << endl;
  double RESresy = (dresy2*ampresy2+dresy3*ampresy3) / (ampresy2+ampresy3);	
  std::cout << "== " << (dresy2*ampresy2+dresy3*ampresy3) / (ampresy2+ampresy3) << std::endl;
  char titleresy2[560];sprintf(titleresy2,"residualV, #sigma = %f #mum",RESresy*10000.);
  hresy->SetTitle(titleresy2);
  hresy->GetXaxis()->SetTitle("Res. (cm)");


  /////////////////////////////////////////////////////////////////////////////////////////
  cxX->cd(3);hresz->Draw();
  TF1* fgpresz2 = new TF1("fgpresz2", "gaus(0)");
  fgpresz2->SetParName(0, "Ag");
  double max24 = hresz->GetMaximum();
  fgpresz2->SetParameter(0, max24);
  fgpresz2->SetParName(1, "mean1");
  fgpresz2->SetParameter(1, 0.);
  fgpresz2->SetParName(2, "sigma1");
  fgpresz2->SetParameter(2, 0.015);
  hresz->Fit("fgpresz2","","",-0.03,0.03);
  hresz->GetXaxis()->SetRangeUser(-0.05,0.05);
  
  TF1* fggpresz2 = new TF1("fggpresz2", "gaus(0) +  gaus(3)");
  fggpresz2->SetParName(0, "Ag");
  fggpresz2->SetParameter(0, fgpresz2->GetParameter(0));
  fggpresz2->SetParName(1, "mean1");
  fggpresz2->SetParameter(1, fgpresz2->GetParameter(1) );
  fggpresz2->SetParName(2, "sigma1");
  fggpresz2->SetParameter(2, fabs(fgpresz2->GetParameter(2)));
  fggpresz2->SetParName(3, "Ag2");
  fggpresz2->SetParameter(3, max24);
  fggpresz2->SetParName(4, "mean2");
  fggpresz2->SetParameter(4, 0.);
  fggpresz2->SetParName(5, "sigma2");
  fggpresz2->SetParameter(5, 0.1);
  fggpresz2->SetLineWidth(2); fggpresz2->SetNpx(1000);
  hresz->GetXaxis()->SetRangeUser(-0.2,0.2);
  hresz->Fit("fggpresz2","","",-0.2,0.2);
  hresz->Draw();
  //draw gaussians
  TF1* fg1resz2 = new TF1("fg1resz2", "gaus", -0.3,0.3);
  for(int i = 0;i<3;++i) fg1resz2->SetParameter(i,fggpresz2->GetParameter(i));
  fg1resz2->SetLineColor(kBlue);fg1resz2->SetLineWidth(2); fg1resz2->SetNpx(1000);
  fg1resz2->Draw("same");
  
  TF1* fg2resz2 = new TF1("fg2resz2","gaus", -0.3,0.3);
  for(int i = 0;i<3;++i) fg2resz2->SetParameter(i,fggpresz2->GetParameter(i+3));
  fg2resz2->SetLineColor(kRed);fg2resz2->SetLineWidth(2); fg2resz2->SetNpx(1000);
  fg2resz2->Draw("same");
 
  double dresz2 = fabs(fggpresz2->GetParameter(2));double ampresz2 = fggpresz2->GetParameter(0);
  double dresz3 = fabs(fggpresz2->GetParameter(5));double ampresz3 = fggpresz2->GetParameter(3);
  std::cout << dresz2 << " x " << ampresz2 << "\t" << dresz3 << " x " << ampresz3 << endl;
  double RESresz = (dresz2*ampresz2+dresz3*ampresz3) / (ampresz2+ampresz3);	
  std::cout << "== " << (dresz2*ampresz2+dresz3*ampresz3) / (ampresz2+ampresz3) << std::endl;
  char titleresz2[560];sprintf(titleresz2,"residualW (3.5<W<4), #sigma = %f #mum",RESresz*10000.);
  hresz->SetTitle(titleresz2);
  hresz->GetXaxis()->SetTitle("Res. (cm)");

}//end of macro


Double_t diffufit(Double_t * x , Double_t*par)
{
  double z = x[0]+par[3];
  double reso =par[0];
  double Dt =par[1];
  double Neff =par[2];
  double miaou = reso*reso + Dt*Dt*z/Neff;
  if(miaou>0) return sqrt(miaou);
  else return 0.;


}