// sample to plot ND data
// 04/26/2007
//
//
// C/C++ headers
//
#include "TrAng.h"
#include "QtoInc.h"
#include "option.h"

#include <fstream>
#include <iostream>
#include <stdio.h>
//
// ROOT headers
//
#include <TTree.h>
#include <TObject.h>
#include <TList.h>
#include <TArrayI.h>
#include <TSystem.h>
#include <TSystemDirectory.h>
#include <TString.h>
#include <TFile.h>
#include <TCanvas.h>
//#include <TMatrixD.h>
#include "TStyle.h"
//#include <TH2F.h>
#include <TGraphErrors.h>
#include <TMultiGraph.h>

#include <PamLevel2.h>
#include <TStyle.h>

using namespace std;

TH2F* NormalizeTH2F(TH2F* h1,Int_t Nev){
    TH2F* h2 = new TH2F("h2","Zenit, deg; Azimut, deg",36,0,360,20,140,180);
    //h2-SetName("h2");
    //h2.SetBins(36,0,360,20,140,180);
    for(Int_t i=0;i<=36;i++) for(Int_t j=0;j<=20;j++)h2->Fill(i*10-5,j*2+139,h1->GetBinContent(i,j)/Nev);
    return h2;
}

void DivideHystTH2F(TString file1, TString file2, TString hyst1, TString hyst2, Int_t Nev1, Int_t Nev2, string outdir,Bool_t Err){
    gStyle->SetPalette(1);
    TFile* f1 = new TFile(file1);
    TFile* f2 = new TFile(file2);
    TH2F* h1 = (TH2F*)f1->Get(hyst1);
    TH2F* h2 = (TH2F*)f2->Get(hyst2);
    TH2F* h1n = NormalizeTH2F(h1,1);
    TH2F* h2n = NormalizeTH2F(h2,1);
    /*TH2F* h1n = new TH2F("h1n","Zenit, deg; Azimut, deg",36,0,360,20,140,180);
    TH2F* h2n = new TH2F("h2n","Zenit, deg; Azimut, deg",36,0,360,20,140,180);
    for(Int_t i=0;i<=36;i++) for(Int_t j=0;j<=20;j++){
	h1n->Fill(i*10-5,j*2+140,h1->GetBinContent(i,j)/h1->GetEntries());
	h2n->Fill(i*10-5,j*2+140,h2->GetBinContent(i,j)/h2->GetEntries());
    }*/
    //vector<TGraphErrors> GrafErr(8);
    TH2F* h3 = new TH2F("h3","Zenit angle, deg; Azimut angle, deg",36,0,360,20,140,180);
    h3->Divide(h1n,h2n);
    TH2F* h3ErrPer;
    if(Err){
    	TH2F* h3Err = new TH2F("h3Err","Errot of Zenit angle, deg; Azimut angle, deg",36,0,360,20,140,180);
        h3ErrPer = new TH2F("h3Err","Errot of Zenit angle, deg; Azimut angle, deg",36,0,360,20,140,180);
    	//ofstream fe((TString)outdir+"Divide_errors"+hyst1+"_"+hyst2+".txt");
    	ofstream fe2((TString)outdir+"Protons_Divide_errors"+hyst1+"_"+hyst2+".txt");
        //Int_t Nev1 = 0;
    	//Int_t Nev2 = 0;
        //for(Int_t i=1;i<=36;i++) for(Int_t j=0;j<=20;j++){Nev1=Nev1+h1->GetBinContent(i,j);Nev2=Nev2+h2->GetBinContent(i,j);}
    	//cout<<"Nev1 = "<<Nev1<<" Nev2 = "<<Nev2<<endl;
        
        for(Int_t j=20;j>=1;j--){
    	    for(Int_t i=1;i<=36;i++){
		//fe<<i<<"\t"<<j<<"\t"<<h1->GetBinContent(i,j)<<"/"<<h2->GetBinContent(i,j)<<" = ";
	        if(h2->GetBinContent(i,j)!=0){
		    Float_t Val = (h1->GetBinContent(i,j))/(h2->GetBinContent(i,j));
		    Float_t Err = 0;
		    if(h1->GetBinContent(i,j)==0)Err=0;else Err = Val*(sqrt(1./h1->GetBinContent(i,j)+1./h2->GetBinContent(i,j)));
		    //Float_t Err = (h2->GetBinContent(i,j)*sqrt(h1->GetBinContent(i,j))-h1->GetBinContent(i,j)*sqrt(h2->GetBinContent(i,j)))/pow(h2->GetBinContent(i,j),2);
		    //Float_t Err = (h1->GetBinContent(i,j)/pow(h2->GetBinContent(i,j),2))+pow(h1->GetBinContent(i,j),2)*h2->GetBinContent(i,j)/pow(h2->GetBinContent(i,j),4);
		    Float_t ErrPer = Err*100/(fabs(Val-1));
		    fe2<<(-j+20)*2+1<<"\t"<<i*10-5<<"\t"<<h1->GetBinContent(i,j)<<"\t"<<h2->GetBinContent(i,j)<<endl;
		    //fe<<Val<<" +/- "<<Err<<"\t";
		    if(ErrPer>100)ErrPer=100;
		    h3ErrPer->Fill(i*10-5,j*2+139,ErrPer);
		    h3Err->Fill(i*10-5,j*2+139,Err);
	        }
	    }
	//fe<<"\n";
	}

	//h1ErrCanv->Divide(1,8);
	for(Int_t i=0;i<8;i++){
	    Float_t x[36];
	    Float_t y[36];
	    Float_t y1[36];
	    Float_t Ex[36];
	    Float_t Ey[36];
	    TH1F* Profil1D = new TH1F("Profil1D","Profil1D",36,0,36);
	    Profil1D->SetFillColor(2);
	    for(Int_t j=0;j<36;j++){
		y[j] = h3->GetBinContent(j,19-i);
	        Ey[j] = h3Err->GetBinContent(j,19-i);
		x[j] = j;
	        Ex[j] = 0;
		y1[j] = 1;
	        Profil1D->Fill(j,h2->GetBinContent(j,19-i));
	    }
	    TCanvas* h1ErrCanv = new TCanvas("h1ErrCanv","",1200,800);
	    h1ErrCanv->Divide(1,2);
	    TGraphErrors* Gep = new TGraphErrors(36,x,y,Ex,Ey);
	    TGraph* _1 = new TGraph(36,x,y1);
	    Gep->SetLineColor(4);
	    Gep->SetFillColor(2);
	    h1ErrCanv->cd(1);
	    Gep->Draw("AB");
	    _1->Draw("");
	    h1ErrCanv->cd(2);
	    Profil1D->Draw("");
	    /*string s;
	    stringstream out;
	    out<<i;
	    s=out.str();
	    h1ErrCanv->SaveAs((TString)outdir+"Divide_"+hyst1+"_"+hyst2+"_1Ds"+s+".png");*/
	    Gep->Delete();
	    delete h1ErrCanv;
	}
	/*fe2.close();
        TFile f3((TString)outdir+"DivideHyst"+hyst1+"_"+hyst2+".root","recreate");
	h3->Write();
        f3.Close();*/
    }    
    /*TCanvas* c1 = new TCanvas("c1",hyst1+"/"+hyst2,1200,1200);
    if(Err) c1->Divide(1,2);
    c1->cd(1);
    h3->Draw("colz");
    if(Err){
	c1->cd(2);
        h3ErrPer->Draw("colz");
    }
    c1->SaveAs((TString)outdir+"Divide_"+hyst1+"_"+hyst2+".png");
    delete c1;*/
    h2->Delete();
    h1->Delete();
    h3->Delete();
    f1->Delete();
    f2->Delete();
}

void DiferenceHystTH2F(TString file1, TString file2, TString hyst1, TString hyst2, Int_t Nev1, Int_t Nev2, string outdir){
    gStyle->SetPalette(1);
    TFile* f1 = new TFile(file1);
    TFile* f2 = new TFile(file2);
    TH2F* h1 = (TH2F*)f1->Get(hyst1);
    TH2F* h2 = (TH2F*)f2->Get(hyst2);
    TH2F* h1n = NormalizeTH2F(h1,Nev1);
    TH2F* h2n = NormalizeTH2F(h2,Nev2);
    TH2F* h3 = new TH2F("h3","Zenit angle, deg; Azimut angle, deg",36,0,360,20,140,180);
    TH2F* h3Err = new TH2F("h3Err","Zenit angle, deg; Azimut angle, deg",36,0,360,20,140,180);
    TH2F* h3ErrPer = new TH2F("h3Err","Errot of Zenit angle, deg; Azimut angle, deg",36,0,360,20,140,180);
    ofstream fe((TString)outdir+"Diference_errors2.txt");
    for(Int_t i=1;i<=36;i++){
	for(Int_t j=1;j<=20;j++){
	    Float_t Val = h1->GetBinContent(i,j)/Nev1-h2->GetBinContent(i,j)/Nev2;
	    Float_t Err = sqrt(pow(h1->GetBinContent(i,j)/Nev1,2)*(1./h1->GetBinContent(i,j)+1./Nev1)+pow(h2->GetBinContent(i,j)/Nev2,2)*(1./h2->GetBinContent(i,j)+1./Nev2));
	    Float_t ErrPer = Err*100/(fabs(Val));
	    h3->Fill(i*10-5,j*2+139,h1->GetBinContent(i,j)/Nev1-h2->GetBinContent(i,j)/Nev2);
	    fe<<i<<"\t"<<j<<"\t"<<h1->GetBinContent(i,j)/Nev1<<" - "<<h2->GetBinContent(i,j)/Nev2<<" = "<<Val;
	    if(h1->GetBinContent(i,j)!=0 && h2->GetBinContent(i,j)!=0){
		fe<<" +/- "<<Err<<" (+/- "<<ErrPer<<" )\n";
		if(ErrPer>100)ErrPer=100;
		h3Err->Fill(i*10-5,j*2+139,Err);
		h3ErrPer->Fill(i*10-5,j*2+139,ErrPer);
	    }else fe<<"\n";
	}
    fe<<"\n";
    }
    for(Int_t i=0;i<8;i++){
	Float_t x[36];
	Float_t y[36];
	Float_t Ex[36];
	Float_t Ey[36];
	TH1F* Profil1D = new TH1F("Profil1D","Profil1D",36,0,36);
	Profil1D->SetFillColor(2);
	for(Int_t j=0;j<36;j++){
	    if(h1n->GetBinContent(j,19-i)!=0)y[j] = h3->GetBinContent(j,19-i)*100/h1n->GetBinContent(j,19-i); else y[j]=0;
	    if(h1n->GetBinContent(j,19-i)!=0) Ey[j] = h3Err->GetBinContent(j,19-i)*100/h1n->GetBinContent(j,19-i); else Ey[j]=0;
	    Profil1D->Fill(j,h2->GetBinContent(j,19-i));
	    //cout<<"h3 = "<<h3->GetBinContent(j,19-i)<<" h1n = "<<h1n->GetBinContent(j,19-i)<<" yn = "<<y[j]<<endl;
	    //Ey[j] = h3ErrPer->GetBinContent(j,19-i);
	    //cout<<Ey[j]<<endl;
	    x[j] = j;
	    Ex[j] = 0;
	}
	TCanvas* h1ErrCanv = new TCanvas("h1ErrCanv","",1200,800);
	h1ErrCanv->Divide(1,2);
	TGraphErrors* Gep = new TGraphErrors(36,x,y,Ex,Ey);
	Gep->SetLineColor(4);
	Gep->SetFillColor(2);
	h1ErrCanv->cd(1);
	Gep->Draw("AB");
	h1ErrCanv->cd(2);
	Profil1D->Draw("");
	string s;
	stringstream out;
	out<<i;
	s=out.str();
	h1ErrCanv->SaveAs((TString)outdir+"Diference_"+hyst1+"_"+hyst2+"_1Ds"+s+".png");
	Gep->Delete();
	Profil1D->Delete();
	delete h1ErrCanv;
    }
    fe.close();
    TFile f3((TString)outdir+"resultDifferenceHyst.root","recreate");
    h3->Write();
    f3.Close();
    TCanvas* c1 = new TCanvas("c1",hyst1+" - "+hyst2,1200,1200);
    c1->Divide(1,2);
    c1->cd(1);
    h3->Draw("colz");
    c1->cd(2);
    h3ErrPer->Draw("colz");
    c1->SaveAs((TString)outdir+"Difference_"+hyst1+"_"+hyst2+".png");
    delete c1;
    h2->Delete();
    h1->Delete();
    h3->Delete();
    f1->Delete();
    f2->Delete();
}

void MergeHystTH2F(TString file1, TString file2, TString hyst1, TString hyst2, string outdir){
    gStyle->SetPalette(1);
    TFile* f1 = new TFile(file1);
    TFile* f2 = new TFile(file2);
    TH2F* h1 = (TH2F*)f1->Get(hyst1);
    TH2F* h2 = (TH2F*)f2->Get(hyst2);
    TH2F* h3 = new TH2F("h3","Zenit angle, deg; Azimut angle, deg",36,0,360,20,140,180);
    for(Int_t i=0;i<=36;i++) for(Int_t j=0;j<=20;j++) h3->Fill(i*10-5,j*2+140,h1->GetBinContent(i,j)+h2->GetBinContent(i,j));
    TFile f3((TString)outdir+hyst1+"_"+hyst2+"_resultMergeHyst.root","recreate");
    h3->Write();
    f3.Close();
    TCanvas* c1 = new TCanvas("c1",hyst1+" + "+hyst2,1200,800);
    c1->cd();
    h3->Draw("colz");
    c1->SaveAs((TString)outdir+"Merge_"+hyst1+"_"+hyst2+".png");
    delete c1;
    h2->Delete();
    h1->Delete();
    h3->Delete();
    f1->Delete();
    f2->Delete();
}


void Convert2Dto1Dhyst(TString file, TString hyst, string outdir){
    TFile* f1 = new TFile(file);
    TH2F* h1 = (TH2F*)f1->Get(hyst);
    //TH1F* h2 = new TH1F("h2","Zenit Efficiency",20,140,180);
    //TH1F* h3 = new TH1F("h3","Azimuth Efficiency",36,0,360);
    //for(Int_t i=0;i<36;i++) for(Int_t j=0;j<20;j++) h3->Fill(i*10-5,h1->GetBinContent(i,j));
    //for(Int_t i=0;i<20;i++) for(Int_t j=0;j<36;j++) h2->Fill(j*2+140,h1->GetBinContent(j,i));
    TH1D* h2 = (TH1D*)h1->ProjectionY(" ",1,20,"");
    TCanvas* c1 = new TCanvas("c1","Zenit Efficiency",1200,800);
    c1->cd();
    h2->Draw();
    c1->SaveAs((TString)outdir+hyst+"Zenit.png");
    delete c1;
    TFile f2((TString)outdir+hyst+"Zenit.root","recreate");
    h2->Write();
    f2.Close();
    TH1D* h3 = (TH1D*)h1->ProjectionX(" ",1,36,"");
    TCanvas* c2 = new TCanvas("c2","Azimuth Efficiency",1200,800);
    c2->cd();
    h3->Draw();
    c2->SaveAs((TString)outdir+hyst+"Azimuth.png");
    delete c2;
    TFile f3((TString)outdir+hyst+"Azimuth.root","recreate");
    h3->Write();
    f3.Close();
    h3->Delete();
    h1->Delete();
    h2->Delete();
    f1->Delete();
}

void ReBinTH2F(TString file, TString hyst, string outdir){
    gStyle->SetPalette(1);
    TFile* f1 = new TFile(file);
    TH2F* h1 = (TH2F*)f1->Get(hyst);
    TH2F* h2 = new TH2F("h2","Zenith, deg; Azimuth, deg",18,0,360,20,140,180);
    for(Int_t i=0;i<=20;i++) for(Int_t j=0;j<=36;j++){h2->Fill(j*10-5,i*2+140,h1->GetBinContent(j,i)+h1->GetBinContent(j+1,i));j++;} 
    TFile f2((TString)outdir+hyst+"Re-Bin.root","recreate");
    h2->Write();
    f2.Close();
    TCanvas* c1 = new TCanvas("c1",hyst,800,600);
    c1->cd();
    h2->Draw("colz");
    c1->SaveAs((TString)outdir+hyst+"Re-Bin.png");
    delete c1;
    h2->Delete();
    h1->Delete();
    f1->Delete();
}

void Modeling(Double_t R, Int_t Nevents, Float_t zmax, string outdir){
    cout<<"Modelling..."<<endl;
    gStyle->SetPalette(1);
    TH2F* h1 = new TH2F("h1","Azimut,deg;Zenit,deg",36,0,360,20,140,180);
    TH2F* h2 = new TH2F("h2","Azimut,deg;Zenit,deg",36,0,360,20,140,180);
    TH2F* h3 = new TH2F("h3","Azimut,deg;Zenit,deg",36,0,360,20,140,180);
    TH2F* h4 = new TH2F("h4","Azimut,deg;Zenit,deg",36,0,360,20,140,180);
    Float_t A;
    Float_t Zenit;
    Float_t Azimuth;
    for(Int_t i=0;i<Nevents;i++){
	A=rand() % 1000; A=A/1000;
	Zenit=TMath::RadToDeg()*acos(pow(A,0.5));
	//Zenit=3;
	Azimuth=rand() % 360000;
	Azimuth=Azimuth/1000;
	h4->Fill(Azimuth,180-Zenit);
	/*for(Int_t nx=0;nx<24;nx++){
	    for(Int_t ny=0;ny<24;ny++){
		h2->Fill(Azimuth,180-Zenit);
		Float_t x = (nx+0.5)*6.75;
		Float_t y = (ny+0.5)*5.5;
		Float_t k = sin(TMath::DegToRad()*Zenit)*cos(TMath::DegToRad()*Azimuth);
		Float_t l = sin(TMath::DegToRad()*Zenit)*sin(TMath::DegToRad()*Azimuth);
		Float_t m = cos(TMath::DegToRad()*Zenit);
		Float_t xl = x+k*zmax/m;
		Float_t yl = y+l*zmax/m;
		Float_t hl = sqrt(pow(x-xl,2)+zmax);
		Float_t za=0; Float_t xa=0;
		Float_t Rk=1000*0.3*R/0.4;
		//cout<<"Zenith = "<<Zenit<<" Azimuth = "<<Azimuth<<endl;
		//cout<<"x = "<<x<<" y = "<<y<<" Rk = "<<Rk<<" xl = "<<xl<<" yl = "<<yl<<endl;
		if(x-xl==0){
		 za==0;
		 xa==x-Rk;
		}else{
		    //Float_t ko = 0;
		    //Float_t koo = 0;
		    if(x-xl<0) {za = zmax-Rk*cos(atan(zmax/fabs(x-xl)));}
		    if(x-xl>0) {za = zmax+Rk*cos(atan(zmax/fabs(x-xl)));}
		    xa = x-Rk*sin(atan(zmax/fabs(x-xl)));
		}
		//cout<<"xa = "<<xa<<" za = "<<za<<" atan = "<<TMath::RadToDeg()*atan(445.1/fabs(x-xl))<<endl;
		Float_t xll;
		if(Rk>za && ((za>0 && xa+Rk<=162 && xa+Rk>=0) || za<=0 )) xll = xa+sqrt(pow(Rk,2)-pow(za,2)); else xll = 10000;
		if(xll<=162 && xll>=0 && yl<=132 && yl>=0) h1->Fill(Azimuth,180-Zenit);
		//cout<<"xll = "<<xll<<endl;
		//if(xl<=162 && xl>=0 && yl<=132 && yl>=0)cin>>za;
	    }
	}*/
    }/*
    TCanvas* c1 = new TCanvas("c1","Model distribution",1200,800);
    h1->Draw("colz");
    c1->SaveAs((TString)outdir+".png");
    TFile f2((TString)outdir+".root","recreate");
    h1->Write();
    f2.Close();
    */
    TCanvas* c3 = new TCanvas("c3","Model Isotropic distribution",1200,800);
    h4->Draw("colz");
    c3->SaveAs((TString)outdir+"_isotrop.png");
    TFile f4((TString)outdir+"_isotrop.root","recreate");
    h4->Write();
    f4.Close();
    /*
    TFile f3((TString)outdir+"Abs.root","recreate");
    h2->Write();
    f3.Close();
    h3->Divide(h2,h1);
    TCanvas* c2 = new TCanvas("c2","Model distribution",1200,800);
    h3->Draw("colz");
    c2->SaveAs((TString)outdir+"Rate.png");
    delete c1;
    delete c2;*/
    delete c3;
    h1->Delete();
    h2->Delete();
    h3->Delete();
    h4->Delete();
}
     
int main(int argc, char* argv[]){

OptionParam opt1;

string outdir = "/home/pamelaprod/malakhov/QtoInc/Picture/";
Bool_t AnglHist = false;
Bool_t DiffHist = false;
Bool_t PamEff = false;
Bool_t DoTr = false;
Bool_t PamAngTime = false;
Bool_t PamExp = false;

/*vector<TH1F> Hyst(1);
Hyst[0].SetName("ElectronLowEnergy");
Hyst[0].SetBins(36,0,360);
Hyst[0].FillRandom("gaus",5000);
TFile f1("/home/pamelaprod/malakhov/QtoInc/TEST/TEST.root","recreate");
Hyst[0].Write();
f1.Close();
//TFile* f1 = new TFile("/home/pamelaprod/malakhov/QtoInc/TEST/TEST.root");
//TH1F* h1 = (TH1F*)f1->Get("");
//cout<<h1->GetEntries()<<endl;
*/

if(argc<2){
    cout<<"You have to insert at least a file to analisys \nUsing --help for more information\n";
    exit(1);
}
//opt1.opt.PE.Proton.Hyst1.resize(opt1.opt.PE.Proton.Hyst1.size()+1);
//opt1.opt.PE.Proton.Hyst1[0].SetBins(36,0,360);
//opt1.opt.PE.Proton.Hyst1[0].SetName("Test2");
//cout<<opt1.opt.PE.Proton.Hyst1.size()<<endl;
//cout<<opt1.opt.PE.Proton.Hyst1[0].GetName()<<endl;

//vector<TH2F>Hyst2DF(1);
//Hyst2DF[0].SetName("Test1");
//Hyst2DF[0].SetBins(36,0,360,36,0,360);
//Hyst2DF.resize(Hyst2DF.size()+1);
//Hyst2DF[1].SetName("Test2");
//Hyst2DF[1].SetBins(36,0,360,36,0,360);
//cout<<Hyst2DF[0].GetName()<<endl;

if(!strcmp(argv[1],"--help")){
    opt1.helprint();
/*    cout<<"\nUsage \n";
    cout<<"\n QtoInc --WorkDir directory_with_level2_files --OutPuth directory_for_output files [Options] \n";
    cout<<"\n --WorkDir		full puth to the Level2 file \n";
    cout<<" --OutPuth 		full puth for putputh files \n";
    cout<<"\n options are:\n";
    cout<<"--help		print this help and exit \n";
    cout<<"-DoTr		Use DoTrack2 procedure for calculating direction of particle flight into Pamela \n";
    cout<<"-PamEff		Calculating Pamela's angular efficiency \n";
    cout<<"-DiffHist		Getting hystogramm for diference between direction getting using DoTrack2 and using axv & ayv values \n";
    cout<<"-AngHist 		Getting hystogramm for various angles(Pitch Angle, Pamela's Main axis angles, etc) [default] \n";
    cout<<"-PamAngTime		Not Useful now \n";
    cout<<"-PamExp		Calculating exposition for each possible direction in Pamela's aperture \n";
    cout<<"\nFor Example: \n";
    cout<<"\n./QtoInc.exe --WorkDir /data01/_3/704_3/ --OutPuth /home/pamelaprod/malakhov/QtoInc/Picture5/ -DiffHist -DoTr \n\n";
*/    exit(1);
}

/*if(!strcmp(argv[1],"--Normalize")){
    if(4 >= argc){
	cout<<"--Normalize  needs argument\n See --help\n";
	exit(1);
    }
    NormalizeTH2F((TString)argv[2],(TString)argv[3],argv[4]);
    exit(0);
}*/

if(!strcmp(argv[1],"--Model")){
    if(5 >= argc){
	cout<<"--Model  needs argument\n See --help\n";
	exit(1);
    }
    Modeling(atof(argv[2]),atoi(argv[3]),atof(argv[4]),argv[5]);
    exit(0);
}

if(!strcmp(argv[1],"--Diference")){
    if(6 >= argc){
	cout<<"--Diference  needs argument\n See --help\n";
	exit(1);
    }
    
    DiferenceHystTH2F((TString)argv[2],(TString)argv[3],(TString)argv[4],(TString)argv[5],atoi(argv[6]),atoi(argv[7]),argv[8]);
    exit(0);
}

if(!strcmp(argv[1],"--Compare")){
    //gStyle->SetPalette(1);
    //TMultiGraph* mg;
    //vector<TGraphErrors*> Gep;
    Float_t Val[20][36];
    Float_t Err[20][36];
    string tmp;
    char g;
    TString out = "";
    for(Int_t ii=2;ii<argc;ii++){
	if(!strcmp(argv[ii],"-out")){
	    out = (TString)argv[ii+1];
	    ii++;
	}else{
	    if(!strcmp(argv[ii],"-in"));
	    ifstream in(argv[ii+1],ios::in);
	    for(Int_t i=20;i>0;i--){
	    for(Int_t j=1;j<=36;j++){
		in>>tmp;in>>tmp;in>>g;
		if(g=='0'){
		    Val[i-1][j-1]=0;
		    Err[i-1][j-1]=0;
		    Bool_t xx=false;
		    Bool_t yy=false;
		    while(g!='='){
			in>>g;
			if(g=='/'){yy=true;continue;}
			if(yy && g=='0') xx=true; else yy=false;
		    }
		    if(!xx){
			in>>tmp;
			in>>tmp;
			in>>tmp;
		    }
		}else{
		    while(g!='/') in>>g;
		    in>>g;
		    if(g=='0'){
			in>>tmp;
			Val[i-1][j-1]=0;
			Err[i-1][j-1]=0;
		    }else{
			while(g!='=') in>>g;
			in>>Val[i-1][j-1];
			in>>tmp;
			in>>Err[i-1][j-1];
		    }
		}
	    }}
	    ii++;
	    //Val.resize(Val.size()+1);
	    //Val[Val.size()-1][1][1]=0;
	}
    }
    ofstream fv(out+"Val.txt");
    ofstream fe(out+"Err.txt");
    for(Int_t i=20;i>0;i--){
    for(Int_t j=1;j<=36;j++){fv<<Val[i-1][j-1]<<"\t";fe<<Err[i-1][j-1]<<"\t";}
        fv<<endl;fe<<endl;
    }
    fv.close();
    fe.close();
}

if(!strcmp(argv[1],"--Merge")){
    if(6 >= argc){
	cout<<"--Merge  needs argument\n See --help\n";
	exit(1);
    }
    gStyle->SetPalette(1);
    //vector<TString> Hysts;
    //vector<TString> Files;
    vector<TH1F*> h1;
    vector<TH2F*> h2;
    TH2F h2sum;
    TH1F h1sum;
    Bool_t H2 = false;
    TString out = "";
    if(!strcmp(argv[2],"2D")) H2=true;
    for(Int_t i=3;i<argc;i++){
	if(!strcmp(argv[i],"-out")){
	    out = (TString)argv[i+1];
	    i++;
	    //break;
	}else{
	    if(H2){
		h2.resize(h2.size()+1);
		cout<<"h2 resize... New size is "<<h2.size()<<endl;
		TFile* f2 = new TFile((TString)argv[i]);
		//cout<<"file is reading"<<endl;
	        //f2->ls();
	        h2[h2.size()-1] = (TH2F*)f2->Get(argv[i+1]);
		//cout<<argv[i+1]<<endl;
		//cout<<h3->GetNbinsX()<<endl;
		//delete h3;
		//delete f1;
		i++;
	    }else{
		h1.resize(h1.size()+1);
	        TFile* f1 = new TFile((TString)argv[i]);
		h1[h1.size()-1] = (TH1F*)f1->Get(argv[i+1]);
	        i++;
	    }
	}
    }
    cout<<"outputh is "<<out<<endl;
    //cout<<"NBinsX = "<<h2[0]->GetNbinsX()<<" Xmin = "<<h2[0]->GetXaxis()->GetXmin()<<" Xmax = "<<h2[0]->GetXaxis()->GetXmax()<<" NbinsY = "<<h2[0]->GetNbinsY()<<" Ymin = "<<h2[0]->GetYaxis()->GetXmin()<<" Ymax = "<<h2[0]->GetXaxis()->GetXmax()<<endl;
    if(H2)h2sum.SetBins(h2[0]->GetNbinsX(),h2[0]->GetXaxis()->GetXmin(),h2[0]->GetXaxis()->GetXmax(),h2[0]->GetNbinsY(),h2[0]->GetYaxis()->GetXmin(),h2[0]->GetYaxis()->GetXmax()); else 
    h1sum.SetBins(h1[0]->GetNbinsX(),h1[0]->GetXaxis()->GetXmin(),h1[0]->GetXaxis()->GetXmax());
    h1sum.SetName("Proton");
    h2sum.SetName("Proton");
    if(H2)for(Int_t i=0;i<h2.size();i++) h2sum=h2sum+(*h2[i]); else for(Int_t i=0;i<h1.size();i++) h1sum=h1sum+(*h1[i]);
    if(H2){
	TFile fsum(out+".root","recreate");
        h2sum.Write();
	fsum.Close();
        TCanvas* c1 = new TCanvas("c1","resultMerge",1200,800);
	c1->cd();
        h2sum.Draw("colz");
	c1->SaveAs(out+".png");
        delete c1;
    }else{
	TFile fsum(out+".root","recreate");
        h1sum.Write();
	fsum.Close();
        TCanvas* c1 = new TCanvas("c1","resultMerge",1200,800);
	c1->cd();
        h1sum.Draw();
	c1->SaveAs(out+".png");
        delete c1;
    }
    //MergeHystTH2F((TString)argv[2],(TString)argv[3],(TString)argv[4],(TString)argv[5],argv[6]);
    exit(0);
}

if(!strcmp(argv[1],"--PlotEfficiency")){
    vector<TGraphErrors> ge;
    TString out = "";
    Int_t t = 1;
    //TGraphErrors tmp;
    for(Int_t i=2;i<argc;i++){
	if(!strcmp(argv[i],"-out")){
	    out = (TString)argv[i+1];
	    i++;
	}
	if(!strcmp(argv[i],"-in")){
	    ifstream inVal((TString)argv[i+1]+"Val.txt",ios::in);
	    ifstream inErr((TString)argv[i+1]+"Err.txt",ios::in);
	    Float_t Val[20][36];
	    Float_t Err[20][36];
	    for(Int_t j=0;j<20;j++)for(Int_t k=0;k<36;k++){inVal>>Val[j][k];inErr>>Err[j][k];}
	    if(!strcmp(argv[i+2],"-yline")){
		Int_t n = atoi(argv[i+3]);
		Float_t x[36];
		Float_t y[36];
		Float_t xe[36];
		Float_t ye[36];
		ge.resize(ge.size()+1);
		//ge[ge.size()-1]->Set(36);
		for(Int_t j=0;j<36;j++){
		    x[j]=j*10+5;
		    xe[j]=5;
		    y[j]=Val[n][j];
		    ye[j]=Err[n][j];
		    ge[ge.size()-1].SetPoint(j,x[j],y[j]);
		    ge[ge.size()-1].SetPointError(j,xe[j],ye[j]);
		    //tmp.SetPoint(j,x[j],y[j]);
		}
		ge[ge.size()-1].SetLineColor(t);
		ge[ge.size()-1].SetTitle((TString)argv[i+1]+" Pitch  "+(TString)argv[i+3]);
		t++;
		//TGraphErrors *tmp = new TGraphErrors(36,x,y,xe,ye);
		
		
		//TGraphErrors df;
		//df = *tmp;
		//cout<<"check"<<endl;
		//ge[ge.size()-1] = tmp;
		//cout<<i<<endl;
		//tmp->Delete();
		i+=3;
	    }
	    if(!strcmp(argv[i+2],"-xline")){
		Int_t p = 0;
		Int_t n=0;
		if(!strcmp(argv[i+3],"All")) p=36;else n = atoi(argv[i+3]);
		Float_t x[20];
		Float_t y[20];
		Float_t xe[20];
		Float_t ye[20];
		for(Int_t j=0;j<20;j++){y[j]=0;ye[j]=0;}
		ge.resize(ge.size()+1);
		for(Int_t j=0;j<20;j++){
		    x[j]=j*2+1;
		    xe[j]=1;
		    for(Int_t yi=0;yi<p;yi++){
			if(p!=0)n=yi;
			y[j]+=Val[j][n];
		        ye[j]+=Err[j][n];
		    }
		    ge[ge.size()-1].SetPoint(j,x[j],y[j]);
		    ge[ge.size()-1].SetPointError(j,xe[j],ye[j]);
		}
		ge[ge.size()-1].SetLineColor(t);
		ge[ge.size()-1].SetTitle((TString)argv[i+1]+" Azimuth  "+(TString)argv[i+3]);
		t++;
		//TGraphErrors* tmp = new TGraphErrors(20,x,y,xe,ye);
		//ge.resize(ge.size()+1);
		//ge[ge.size()-1] = (*tmp);
		//tmp->Delete();
		i+=3;
	    }
	}
    }
    TMultiGraph* mg = new TMultiGraph();
    cout<<"check"<<endl;
    for(Int_t i=0;i<ge.size();i++){
	//ge[i]->SetLineColor(i+1);
	//cout<<i<<endl;
	mg->Add(&ge[i]);
    }
    TCanvas* c1 = new TCanvas("c1","eryt",1200,800);
    c1->cd();
    mg->Draw("AP");
    cout<<"check"<<endl;
    cout<<"out = "<<out<<endl;
    c1->SaveAs(out+".png");
    delete c1;
    mg->Delete();
    exit(1);
}

if(!strcmp(argv[1],"--Divide")){
    if(argc<7){
	cout<<"--Divide  needs argument\n See --help\n";
	exit(1);
    }/*
    Bool_t H2 = false;
    TString out = "";
    if(!strcmp(argv[2],"2D")) H2=true;
    out = (TString)argv[7];
    gStyle->SetPalette(1);
    TFile* f1 = new TFile((TString)argv[3]);
    TFile* f2 = new TFile((TString)argv[5]);
    TH1F* h11D; TH1F* h21D;
    TH2F* h12D; TH2F* h22D;
    TH1F h1n1D; TH1F h2n1D;
    TH2F h1n2D; TH2F h2n2D;
    TH1F h1div; TH2F h2div;
    cout<<"out is "<<out<<endl;
    if(H2){
	h12D = (TH2F*)f1->Get(argv[4]);
        h22D = (TH2F*)f2->Get(argv[6]);
	h1n2D.SetBins(h12D->GetNbinsX(),h12D->GetXaxis()->GetXmin(),h12D->GetXaxis()->GetXmax(),h12D->GetNbinsY(),h12D->GetYaxis()->GetXmin(),h12D->GetYaxis()->GetXmax());
	h2n2D.SetBins(h12D->GetNbinsX(),h12D->GetXaxis()->GetXmin(),h12D->GetXaxis()->GetXmax(),h12D->GetNbinsY(),h12D->GetYaxis()->GetXmin(),h12D->GetYaxis()->GetXmax());
	h2div.SetBins(h12D->GetNbinsX(),h12D->GetXaxis()->GetXmin(),h12D->GetXaxis()->GetXmax(),h12D->GetNbinsY(),h12D->GetYaxis()->GetXmin(),h12D->GetYaxis()->GetXmax());
	for(Int_t i = 0; i < h12D->GetNbinsX(); i++){
	    for(Int_t j = 0; j < h12D->GetNbinsY(); j++){
		h1n2D.SetBinContent(i,j,h12D->GetBinContent(i,j)/h12D->GetEntries());
		h2n2D.SetBinContent(i,j,h22D->GetBinContent(i,j)/h22D->GetEntries());
	    }
	}
        h2div.Divide(&h1n2D,&h2n2D);
    }else{
	h11D = (TH1F*)f1->Get((TString)argv[4]);
        h21D = (TH1F*)f2->Get((TString)argv[6]);
        h1n1D.SetBins(h11D->GetNbinsX(),h11D->GetXaxis()->GetXmin(),h11D->GetXaxis()->GetXmax());
        h2n1D.SetBins(h21D->GetNbinsX(),h21D->GetXaxis()->GetXmin(),h21D->GetXaxis()->GetXmax());
        h1div.SetBins(h21D->GetNbinsX(),h21D->GetXaxis()->GetXmin(),h21D->GetXaxis()->GetXmax());
        for(Int_t i = 0; i<h11D->GetNbinsX();i++){
    	    h1n1D.SetBinContent(i,h11D->GetBinContent(i)/h11D->GetEntries());
    	    h2n1D.SetBinContent(i,h21D->GetBinContent(i)/h21D->GetEntries());
        }
        h1div.Divide(&h1n1D,&h2n1D);
    }
    h2div.SetName("Proton");
    h1div.SetName("Proton");
    if(H2){
	TFile fdiv(out+".root","recreate");
        h2div.Write();
	fdiv.Close();
        TCanvas* c1 = new TCanvas("c1","resultDivide",1200,800);
	c1->cd();
        h2div.Draw("colz");
	c1->SaveAs(out+".png");
        delete c1;
    }else{
	TFile fdiv(out+".root","recreate");
        h1div.Write();
	fdiv.Close();
        TCanvas* c1 = new TCanvas("c1","resultDivide",1200,800);
	c1->cd();
        h1div.Draw();
	c1->SaveAs(out+".png");
        delete c1;
    }*/
    Bool_t Err = false;
    if((argc==10) && !strcmp(argv[9],"-Err")) Err = true;
    DivideHystTH2F((TString)argv[2],(TString)argv[3],(TString)argv[4],(TString)argv[5],atoi(argv[6]),atoi(argv[7]),argv[8],Err);
    exit(0);
}

if(!strcmp(argv[1],"--ReBin")){
    if(4 >= argc){
	cout<<"--ReBin  needs argument\n See --help\n";
	exit(1);
    }
    ReBinTH2F((TString)argv[2],(TString)argv[3],argv[4]);
    exit(0);
}

if(!strcmp(argv[1],"--2Dto1D")){
    if(4 >= argc){
	cout<<"--2Dto1D  needs argument\n See --help\n";
	exit(1);
    }
    Convert2Dto1Dhyst((TString)argv[2],(TString)argv[3],argv[4]);
    exit(0);
}

if(!strcmp(argv[3],"--OutPuth")){
    if (4 >= argc){
    cout<<"--OutPath needs argument\n See --help\n";    
    }else{
	outdir = argv[4];
	opt1.outdir = argv[4];
    }
}else{
    cout<<"\n ERROR:\tYou need --OutPuth as second argument for run programm\n See --help for more information\n"<<endl;
    exit(1);
}

if(!strcmp(argv[1],"--WorkDir")){
    if (2 >= argc){
	cout<<"--WorkDir needs argument\n See --help\n";
	exit(1);
    }
}else{
    cout<<"\n ERROR:\tYou need --WorkDir as first argument for run programm\n See --help for more information\n"<<endl;
    exit(1);
}

if(argc==5){
    //opt1.Default();
    cout<<"MESSAGE:\tSetting default options \n\n";
}

for(Int_t i=5; i<argc; i++){

    if(!strcmp(argv[i],"-do")) opt1.debug = true;
    //opt1.AnglHystFunc(argc,argv);
    //exit(1);
    /*
    Int_t NHyst1D = 0; Int_t NHyst2D = 0;
    if(!strcmp(argv[i],"--AnglHyst")){
	
	if(opt1.debug) cout<<"\n MESSAGE:\t AnglHyst was inicialized \n\n";
	if((i+1)>=argc){
	    //opt1.AnglHystDefault(); 
	    if(opt1.debug) cout<<"Setting default value for AnglHyst option.\n\n";
	}else{
	    
	    Int_t i0 = i+1;
	    
	    Bool_t inside = false; //Check that we go inside loop
	    
	    Int_t NPth = 0; Int_t NMA = 0; Int_t NPNO = 0;
	    opt1.AnglHystFunc(argc,argv);
	    
	}	//if((i+1)>=argc) {...}; else

    }	//AnglHyst
    
    if(!strcmp(argv[i],"--DiffHyst")){
	//opt1.opt.DH.DiffHystOn = true;
	//opt1.DiffHystFunc(argc,argv);
    }*/
    if(!strcmp(argv[i],"--PamEff")){
	//opt1.opt.PE.PamEffOn = true;
	opt1.PamEffFunc(argc,argv);
	/*if(opt1.debug){
	    for(Int_t i = 0; i<opt1.Hyst2D.size();i++){
		cout<<"\t\tName of Hystogramm - "<<opt1.Hyst2D[i].Hyst2DF.GetName()<<endl;
		cout<<"\t\tRoot = "<<opt1.Hyst2D[i].HI.root<<endl;
		cout<<"\t\tpng = "<<opt1.Hyst2D[i].HI.png<<endl;
		cout<<"\t\tExpType = "<<opt1.Hyst2D[i].HI.ExpType<<endl;
		cout<<"\t\tFile name = "<<opt1.Hyst2D[i].HI.filename<<endl;
	    }
	}*/
    }	//PamEff
    
    if(!strcmp(argv[i],"-DoTr")) opt1.opt.DoTr = true;
    if(!strcmp(argv[i],"-v") || !strcmp(argv[i],"-verbose")){
	if(i+1>=argc || !atoi(argv[i+1])){
	    cout<<"\n ERROR! \t -verbose needs arguments \n See --help \n";
	    if(opt1.debug) cout<<"Argument number\t"<<i<<"\n";
	    exit(1);
	}
	opt1.opt.verbose = atoi(argv[i+1]);
	if(opt1.opt.verbose>3){
	    cout<<"\n ERROR! \t -verbose needs 0, 1, 2 or 3 only \n See --help \n";
	    if(opt1.debug) cout<<"Argument\t"<<opt1.opt.verbose<<"\n";
	    exit(1);
	}
    }
    if(!strcmp(argv[i],"-s") || !strcmp(argv[i],"-silent")) opt1.opt.verbose = 0;
/*    if(!strcmp(argv[i],"-outform") || !strcmp(argv[i],"-of")){
	Int_t i7 = i+1;
	if(opt1.debug) cout<<"OutForm option was selected...\n\n";
	if(i7>=argc || (strcmp(argv[i7],"jpg") && strcmp(argv[i7],"png") && strcmp(argv[i7],"bmp") && strcmp(argv[i7],"root") && 
	strcmp(argv[i7],"txt") && strcmp(argv[i7],"all"))){
	    cout<<"\n ERROR! \t -outform needs argument \n See --help \n";
	    if(opt1.debug) cout<<"Argument number - "<<i7<<" argc = "<<argc<<"\n";
	    exit(1);
	}
	while(i7<argc && (!strcmp(argv[i7],"jpg") || !strcmp(argv[i7],"png") || !strcmp(argv[i7],"bmp") || !strcmp(argv[i7],"root") || 
	!strcmp(argv[i7],"txt") || !strcmp(argv[i7],"all"))){
	    if(!strcmp(argv[i7],"jpg")) opt1.opt.ext[0] = true;
	    if(!strcmp(argv[i7],"png")) opt1.opt.ext[1] = true;
	    if(!strcmp(argv[i7],"bmp")) opt1.opt.ext[2] = true;
	    if(!strcmp(argv[i7],"root")) opt1.opt.ext[3] = true;
	    if(!strcmp(argv[i7],"txt")) opt1.opt.ext[4] = true;
	    if(!strcmp(argv[i7],"all")){
		opt1.opt.ext[1] = true;
		opt1.opt.ext[3] = true;
		opt1.opt.ext[4] = true;
	    }
	i7++;
	}
    }*/
    if(!strcmp(argv[i],"-ne") || (!strcmp(argv[i],"-nevents"))){
	if(i+1>=argc || !atoi(argv[i+1])){
	    cout<<"\n ERROR! \t -nevents  needs ULong integer \n See --help \n";
	    exit(1);
	}
	if(atoi(argv[i+1])<0){
	    cout<<"\n ERROR! \t -nevents  needs ULong integer \n See --help \n";
	    exit(1);
	}
	opt1.opt.Neve = atoi(argv[i+1]);
    }
    if(!strcmp(argv[i],"-nv") || (!strcmp(argv[i],"-Nverbose"))){
	if(i+1>=argc || !atoi(argv[i+1])){
	    cout<<"\n ERROR! \t -nevents  needs ULong integer \n See --help \n";
	    exit(1);
	}
	if(atoi(argv[i+1])<0){
	    cout<<"\n ERROR! \t -nevents  needs ULong integer \n See --help \n";
	    exit(1);
	}
	opt1.opt.Nverb = atoi(argv[i+1]);
    }
    if(!strcmp(argv[i],"-ns") || (!strcmp(argv[i],"-Nstart"))){
	if(i+1>=argc || !atoi(argv[i+1])){
	    cout<<"\n ERROR! \t -Nstart needs ULong integer \n See --help \n";
	    exit(1);
	}
	if(atoi(argv[i+1])<0){
	    cout<<"\n ERROR! \t -nstart  needs ULong integer \n See --help \n";
	    exit(1);
	}
	opt1.opt.Nstart = atoi(argv[i+1]);
    }
    if(!strcmp(argv[i],"--Show")){
	Int_t a=0;
	if(opt1.debug) a=1;
	if(i>5+a){
	    cout<<"\n Error!: \t --Show must be used as single option without any other \n See --help \n";
	    if(opt1.debug) cout<<"Argument number - "<<i<<" argc = "<<argc<<" arg0 - "<<argv[0]<<" arg1 - "<<argv[1]<<"\n";
	    exit(1);
	}
	if(opt1.debug) cout<<"\nShow Option was selected...\n\n";
	opt1.opt.SP.Existed = true;
	if(!strcmp(argv[i+1],"-nevents")) opt1.opt.SP.Nev = true;
	if(!strcmp(argv[i+1],"-StartTime")) opt1.opt.SP.ST = true;
	if(!strcmp(argv[i+1],"-EndTime")) opt1.opt.SP.ET = true;
	if(!strcmp(argv[i+1],"-FullInformation")){
	    Int_t j6 = i+2;
	    if(opt1.debug) cout<<"\nFull Information option was selected...\n\n";
	    if(j6>=argc){
		cout<<"\n ERROR! \t -FullInformation needs argument \n See --help \n";
		exit(1);
	    }
	    while(j6<argc && (!strcmp(argv[j6],"-absTime") || !strcmp(argv[j6],"-event") || !strcmp(argv[j6],"-UTS") || !strcmp(argv[j6],"-Track"))){
		if(!strcmp(argv[j6],"-event")){
		    if(j6+1>=argc || !atoi(argv[j6+1])){
			cout<<"\n ERROR! \t -nevents needs argument \n See --help \n";
			exit(1);
		    }
		    opt1.opt.SP.FI.NEvent = atoi(argv[j6+1]);
		    if(opt1.debug) cout<<"Event number\t"<<opt1.opt.SP.FI.NEvent<<"\n";
		    j6++;
		}
		if(!strcmp(argv[j6],"-absTime")){
		    if(j6+1>=argc || !atoi(argv[j6+1])){
			cout<<"\n ERROR! \t -absTime needs argument \n See --help \n";
			exit(1);
		    }
		    opt1.opt.SP.FI.AbsTime = atoi(argv[j6+1]);
		    if(opt1.debug) cout<<"Absolute time of Event\t"<<opt1.opt.SP.FI.AbsTime<<"\n";
		    j6++;
		}
		if(!strcmp(argv[j6],"-UTS")){
		    if(j6+1>=argc){
			cout<<"\n ERROR! \t -UTS needs argument \n See --help \n";
			exit(1);
		    }
		    opt1.opt.SP.FI.UTS = argv[j6+1];
		    if(opt1.debug) cout<<"UTS of Event\t"<<opt1.opt.SP.FI.NEvent<<"\n";
		    j6++;
		}
		if(!strcmp(argv[j6],"-Track")) opt1.opt.SP.FI.WT = true;
	    j6++;
	    }
	}
	if(!strcmp(argv[i+1],"-absTime")){
	    if(i+2>=argc || !atoi(argv[i+2])){
		cout<<"\n ERROR! \t -event needs integer argument \n See --help \n";
		exit(1);
	    }
	    opt1.opt.SP.AbsTime = atoi(argv[i+2]);
	    if(opt1.debug) cout<<"AbsTime for convert is "<<opt1.opt.SP.AbsTime<<"\n";
	}
	if(!strcmp(argv[i+1],"-UTS")){
	    if(i+2>=argc){
		cout<<"\n ERROR! \t -UTS needs argument \n See --help \n";
		exit(1);
	    }
	    opt1.opt.SP.UTS = argv[i+2];
	    if(opt1.debug) cout<<"UTS for convert is "<<opt1.opt.SP.UTS<<"\n";
	}
    }
}

//Bool_t OptDef = false;
//for(Int_t t=0; t<5; t++) if (opt1.opt.ext[t] == true) OptDef = true;
//if(!OptDef) opt1.opt.ext[1] = true;
/*
if(opt1.debug){
    if(opt1.opt.SP.Existed){
	cout<<"\nShow options are:\n\n";
	cout<<"number of events - "<<opt1.opt.SP.Nev<<";\t Start time - "<<opt1.opt.SP.ST<<";\t End time - "<<opt1.opt.SP.ET<<";\n";
    }else{
        cout<<"\nOptions are:\n\n";
	cout<<"DoTrack - "<<opt1.opt.DoTr<<";\t verbose - "<<opt1.opt.verbose<<";\n";
	cout<<"output file:\t";
	for(Int_t t=0; t<5; t++) cout<<opt1.opt.ext[t]<<",  ";
	cout<<"\n\n";
    }
}*/

if(!opt1.debug) Calculate(argv[2],outdir,AnglHist,DiffHist,PamEff,DoTr,PamAngTime,PamExp,opt1);

}

//Function to Convert Rigidity to Energy

Float_t ConvertR2T(Float_t R, Float_t M, Float_t Z)
{
 //
 // Convert rigidity (GV) in kinetic energy (GeV) per nucleon
 // input = rigidity (GV), mass (Gev/c^2), A,Z.
 //

 Float_t EcinperN = (sqrt(pow((Z*R),2)+pow(M,2))-M);
 return EcinperN;
}

//Function to convert Energy to Rigidity

Float_t ConvertT2R(Float_t T, Float_t M, Float_t A, Float_t Z)
{
 //
 // Convert kinetic energy (GeV) per nucleon in rigidity (GV)  
 // output = rigidity (GV),input kin energy mass (Gev/c^2), A,Z.
 //

 Float_t R= (1/Z)*(sqrt(pow((A*T+M),2)-pow(M,2)));
 return R;
}


void Calculate(char* dirname, string outdir, Bool_t AnglHist, Bool_t DiffHist, Bool_t PamEff, Bool_t DoTr, Bool_t PamAngTime, Bool_t PamExp, OptionParam opt){

    /**********************************************/
    //First initialization, general for all purposes
    /**********************************************/
    
    if(opt.opt.verbose>=2) cout<<"\nGeneral variables initialisation...\n\n";
    
    TSystemDirectory* workdir = new TSystemDirectory("workdir",dirname);
    //TSystemDirectory* workdir2 = new TSystemDirectory("workdir2","/data01/_3/607_3/");
    TList* flist = workdir->GetListOfFiles();
    PamLevel2* pam_events = new PamLevel2();
    PamelaOrientation* PO = new PamelaOrientation();
    TTree *T = pam_events->GetPamTree(flist,"treename");
    ULong_t nevents = T->GetEntries();

    if(opt.opt.verbose>=1) cout<<"\n Nevents = "<<nevents<<"\n\n";
    
    
    /********************************************************************************/
    /*****NEED TO  CHANGE FOR OTHER COMPUTERS****************************************/
    pam_events->GetTrkLevel2()->LoadField("/data01/pam7/installed/pamela_software/calib/trk-param/field_param-0/");
    /********************************************************************************/
    
    Int_t nz = 6; Float_t zin[6];
    for (Int_t ip=0;ip<nz;ip++) {zin[ip] = pam_events->GetToFLevel2()->GetZTOF(pam_events->GetToFLevel2()->GetToFPlaneID(ip));cout<<zin[ip]<<endl;}
    Trajectory *tr = new Trajectory(nz,zin);

    if(opt.opt.verbose>=2) cout<<"\n Other Variables initialisation...\n\n";

    PamTrack *track;
    
    Int_t ntr;			Float_t Argv = 0;	Double_t PamAzim = 0;
    Double_t PamZenit = 0;	Double_t MyAzim = 0;	Double_t MyZenit = 0;
    Double_t Px = 0;		Double_t Py = 0;	Double_t Pz = 0;
    
    //Peremennye dlya opredeleniya pervogo vhoda v cycly v PamExp
    
    ULong_t PHE1 = 0;		ULong_t PLE1 = 0;	ULong_t PAE1 = 0;
    ULong_t EHE1 = 0;		ULong_t ELE1 = 0;	ULong_t EAE1 = 0;
    ULong_t pHE1 = 0;		ULong_t pLE1 = 0;	ULong_t pAE1 = 0;
    ULong_t AHE1 = 0;		ULong_t ALE1 = 0;	ULong_t AAE1 = 0;
    
    /*********************************************************************************/
    //Histogramms for Count and exposition of Angles (Pitch, Pamela's Main axis, etc )
    /*********************************************************************************/
    
    if(opt.opt.verbose>=2) cout<<"\n--AnglHyst Hystogramm block...\n\n";
    
    /***************************************************************************************/
    //Histogramms for differences in results between calculating using DoTrack2
    //and using axv[0], ayv[0]
    /**************************************************************************************/
    
    if(opt.opt.verbose>=2) cout<<"\n--DiffHyst Hystogramm block...\n\n";

    
    /***************************************************************************************/
    //Initialization histogramms for count of protons depends from Energy. Volodia.
    /***************************************************************************************/
    
//    TH1F* proton10log= new TH1F("Protons Flux","",80,-1.,3.);
//    proton10log->GetXaxis()->SetTitle("log10(Energy) , GeV");
//    TH1F* proton1log= new TH1F("Protons Flux","",80,-1.,3.);
    //TH1F* binw= new TH1F("w","",80,-1.,3.);
    TH1F binw;
    binw.SetBins(80,-1.,3.);
//    proton1log->GetXaxis()->SetTitle("log10(Energy) , GeV");
//    proton1log->SetLineColor(kRed);
//    TCanvas *plogCanvas = new TCanvas("proton flux","protonflux", 800, 800);
    // this is bins wide to calculate Flux
    //Float_t x;
    //Float_t xmin, xmax;
    //Float_t binwide[80];
    for (Int_t l=0 ; l<80; l++) {
	binw.Fill(-1.+0.05*l,pow(10.,0.05*(Float_t)(l+1)-1.)- pow(10.,0.05*(Float_t)(l)-1.)); 
	//binwide[l]=-1.+0.05*l,pow(10.,0.05*(Float_t)(l+1)-1.)- pow(10.,0.05*(Float_t)(l)-1.)
    }
   
    //for (Int_t l=0 ; l<80; l++) {
    //binwide[l]=binw->GetBinContent(l+1);
//    }

    /****************************************************************************/
    //Geometrical Factor. Volodia
    /****************************************************************************/
/*
    TH1F* Geomfactor = new TH1F("G","G",1000,0.1,500.1);
    TH1F* Geomfactorlog= new TH1F("Glog","Glog",80,-1.,3.);
    TH1F* Geomfactorlogelectron= new TH1F("Gloge","Gloge",80,-1.,3.);
    for (Int_t l=0 ; l <80 ;l++) {
	x=pow(10.,(Float_t) 0.05*l-1.);
	Geomfactorlogelectron->SetBinContent(l+1,pow(1/pow(66.7051+50.05404*log10(x),3.5083)+1./pow(21.6238,3.5083),(-1./3.5083)));
	x=ConvertT2R(x,0.938,1., 1.);
	Geomfactorlog->SetBinContent(l+1,pow(1/pow(66.7051+50.05404*log10(x),3.5083)+1./pow(21.6238,3.5083),(-1./3.5083)));
    }//geomfactor for linear scale 
    for (Int_t l=0 ; l <1000 ;l++) {
	x=(Float_t) 0.5*l+0.1;
	Geomfactor->SetBinContent(l+1,pow(1/pow(66.7051+50.05404*log10(x),3.5083)+1./pow(21.6238,3.5083),(-1./3.5083)));
    }
*/
    /****************************************************************************/
    //General loop
    /****************************************************************************/

    if(opt.opt.Neve == 0) opt.opt.Neve = nevents;

    cout<<opt.opt.Neve<<endl;

    for(ULong_t i=opt.opt.Nstart; i<opt.opt.Neve; i++){
    
	//cout<<i<<endl;
	T->GetEntry(i);
	Bool_t OneTrack = false;

	if(pam_events->GetTrkLevel2()->GetNTracks()==1) OneTrack = true;
	
	Int_t M0DE = pam_events->GetOrbitalInfo()->mode;  //0 is zero here
	Float_t Bx = -pam_events->GetOrbitalInfo()->Bdown;			//don't need for PamExp ExpOnly for all geography areas
	Float_t By = pam_events->GetOrbitalInfo()->Beast;			//don't need for PamExp ExpOnly for all geography areas
	Float_t Bz = pam_events->GetOrbitalInfo()->Bnorth;			//don't need for PamExp ExpOnly for all geography areas
	Float_t Babs = pam_events->GetOrbitalInfo()->Babs;			//don't need for PamExp ExpOnly for all geography areas
	Float_t L = pam_events->GetOrbitalInfo()->L;				//don't need for PamExp ExpOnly for all geography areas
	Float_t Alt = pam_events->GetOrbitalInfo()->alt;
	Float_t Lat = pam_events->GetOrbitalInfo()->lat;
	Float_t Lon = pam_events->GetOrbitalInfo()->lon;
	Float_t LieveTime = 0.16*pam_events->GetTrigLevel2()->dltime[0];
	Float_t DeathTime = 0.01*pam_events->GetTrigLevel2()->dltime[1];
	Float_t rigev = 0;
	Float_t detr = 0;
	Int_t hm0; Int_t hm1;
	Int_t hs0; Int_t hs1;

	ntr = 0;    

	Double_t A1; Double_t A2; Double_t A3;
	Double_t SB = 1000; Double_t SA; Double_t SC; Double_t ZC;

	if(OneTrack){
	    track = pam_events->GetTrack(0);				//don't need for PamExp option
	    rigev = 1/track->GetTrkTrack()->GetDeflection();		//dont need for AnglHyst, DiffHyst and PamExp options;
	    detr = track->GetTrkTrack()->GetDEDX();
	    hm0 = pam_events->GetAcLevel2()->hitmap[0];
	    hm1 = pam_events->GetAcLevel2()->hitmap[1];
	    hs0 = pam_events->GetAcLevel2()->hitstatus[0];
	    hs1 = pam_events->GetAcLevel2()->hitstatus[1];
	    
	    if(!OneTrack) ntr = -1; else{
		//cout<<"OneTrack"<<endl;
		Float_t ichi2lim = fabs(track->GetTrkTrack()->GetDeflection())*1.85+3.0;
		if ((track->GetTrkTrack()->chi2<=0) || (track->GetTrkTrack()->chi2>=pow(ichi2lim,4.))) ntr=-1;
	        /*for(Int_t ii=1;ii<=12;ii++){
		    if ((track->GetToFTrack()->beta[ii]<0) || (track->GetToFTrack()->beta[ii]>99)) ntr=-1;
		}*/
		if ((track->GetToFTrack()->beta[12]<0) || (track->GetToFTrack()->beta[12]>99)) ntr=-1;
		if((track->GetTrkTrack()->GetNX()<=4)&&(track->GetTrkTrack()->GetNY()<4)) ntr=-1;
		if(!track->GetTrkTrack()->IsInsideCavity()) ntr=-1;
		if(sqrt(pow(Bx,2)+pow(By,2)+pow(Bz,2))<0.001) ntr=-1;
	        if((M0DE!=0)&&(M0DE!=1)&&(M0DE!=6)&&(M0DE!=2)&&(M0DE!=3)&&(M0DE!=8)&&(M0DE!=4)) ntr=-1;
	    }
	    
	    if(ntr!=-1){
	    
	    if(!opt.opt.DoTr){
	    	
	        Double_t Aaxv = TMath::Abs(track->GetTrkTrack()->axv[0])*TMath::DegToRad();
	        Double_t Aayv = TMath::Abs(track->GetTrkTrack()->ayv[0])*TMath::DegToRad();
	        PamZenit = TMath::RadToDeg()*asin(sqrt(pow(sin(Aayv), 2) + pow(sin(Aaxv), 2)));
	    	
	    	Double_t axv = -track->GetTrkTrack()->axv[0] * TMath::DegToRad();
	    	Double_t ayv = -track->GetTrkTrack()->ayv[0] * TMath::DegToRad();
	    	Double_t angle = atan(sin(TMath::Abs(ayv))/sin(TMath::Abs(axv))) * TMath::RadToDeg();
	    	
	    	PamAzim =  360. - angle;
	    	if(axv>=0 && ayv >=0) PamAzim = angle;
	    	if(axv<0 && ayv >0) PamAzim = 180. - angle;
	    	if(axv<0 && ayv <0) PamAzim = 180. + angle;
	    	
	    	PamAzim = PamAzim * TMath::DegToRad();
	    	PamZenit = (180 - PamZenit) * TMath::DegToRad();
	    	Px = sin(axv);
	    	Py = sin(ayv);
	    	Pz = cos(PamZenit);
	    	    
	    	//cout<<"PamAzimuth = "<<PamAzim<<" PamZenit = "<<PamZenit<<endl;
		
	    }
	    
	    /*****************************************************************************/
	    //Calculate Zenit and Azimutal angle and vector of particle flight in Pamela
	    //using DoTrack2 procedure
	    /*****************************************************************************/
	    
	    if(opt.opt.DoTr){
	    
		track->GetTrkTrack()->DoTrack2(tr);
	        Double_t E11x = tr->x[0];
	        Double_t E11y = tr->y[0];
	        Double_t E11z = zin[0];
	        Double_t E22x = tr->x[3];
	        Double_t E22y = tr->y[3];
	        Double_t E22z = zin[3];
	        Double_t norm = sqrt(pow(E22x-E11x,2)+pow(E22y-E11y,2)+pow(E22z-E11z,2));
	        MyAzim = TMath::RadToDeg()*atan(TMath::Abs(E22y-E11y)/TMath::Abs(E22x-E11x));
	        MyAzim = MyAzim;//-180;
	        MyZenit = TMath::RadToDeg()*acos((E22z-E11z)/norm);
	        if(E22x-E11x>=0 && E22y-E11y <0) MyAzim =  360. - MyAzim;
	        if(E22x-E11x>=0 && E22y-E11y >=0) MyAzim = MyAzim;
	        if(E22x-E11x<0 && E22y-E11y >0) MyAzim = 180. - MyAzim;
	        if(E22x-E11x<0 && E22y-E11y <0) MyAzim = 180. + MyAzim;
	        Px = (E22x-E11x)/norm;
	        Py = (E22y-E11y)/norm;
	        Pz = (E22z-E11z)/norm;
	        
	    }
	    	
	    TMatrixD Fij = PO->ECItoGreenwich(PO->QuatoECI(pam_events->GetOrbitalInfo()->q0,pam_events->GetOrbitalInfo()->q1,pam_events->GetOrbitalInfo()->q2,pam_events->GetOrbitalInfo()->q3),pam_events->GetOrbitalInfo()->absTime);
	    TMatrixD Dij = PO->GreenwichtoGEO(pam_events->GetOrbitalInfo()->lat,pam_events->GetOrbitalInfo()->lon,Fij);
	    TMatrixD Iij = PO->ColPermutation(Dij);
		
	    TMatrixD Eij = PO->PamelatoGEO(Iij,Px,Py,Pz);		//Don't necessary for PamExp
	    	
	    A1 = Iij(0,2);
	    A2 = Iij(1,2);
	    A3 = Iij(2,2);
	    	
	    SB = PO->GetPitchAngle(Eij(0,0),Eij(1,0),Eij(2,0),Bx,By,Bz); 	// Pitch angle
	    SA = PO->GetPitchAngle(A1,A2,A3,Bx,By,Bz);			// Angle between Pamela's main axiz and B
	    SC = PO->GetPitchAngle(1,0,0,Bx,By,Bz);				// Angle between direction to zenit and B
	    
	    }
	}
	
	//cout<<"OneTrack"<<endl;
	if(sqrt(pow(Bx,2)+pow(By,2)+pow(Bz,2))<0.001) ntr=-1;else{
	    //cout<<"US.size = "<<opt.US.size()<<endl;
	    for(Int_t jz = 0; jz<opt.PS.size();jz++){
		for(Int_t jy = 0; jy<opt.US.size();jy++){
	    	    for(Int_t jx = 0; jx<opt.US[jy].USP.size();jx++){
	    		//if(SB<400) cout<<"Pitch = "<<SB<<" Babs = "<<Babs<<" L = "<<L<<" Alt = "<<Alt<<endl;
			//cout<<"Fill Hyst1D["<<thg<<"]"<<endl;
	    		if(Babs<opt.US[jy].USP[jx].Bless && Babs>opt.US[jy].USP[jx].Bmore && L<opt.US[jy].USP[jx].Lless && L>opt.US[jy].USP[jx].Lmore && Alt<opt.US[jy].USP[jx].altless && Alt>opt.US[jy].USP[jx].altmore  && SB<opt.PS[jz].PitchMax && SB>opt.PS[jz].PitchMin){
			    for(Int_t hg = 0; hg<opt.US[jy].NAH1D.size(); hg++){
				Int_t thg = opt.US[jy].NAH1D[hg];
			        for(Int_t ps = 0; ps<opt.PS[jz].NAH1D.size();ps++){
			    	    if(thg==opt.PS[jz].NAH1D[ps]){
			    		opt.Hyst1D[thg].HI.LieveTime+= 0.16*pam_events->GetTrigLevel2()->dltime[0];
				        opt.Hyst1D[thg].HI.DeathTime+= 0.01*pam_events->GetTrigLevel2()->dltime[1];
				        //cout<<"Pitch = "<<SB<<" Babs = "<<Babs<<" L = "<<L<<" Alt = "<<Alt<<endl;
				        //cout<<"Fill Hyst1D["<<thg<<"]"<<endl;
				    }
				}
		    	    }
			    for(Int_t hg = 0; hg<opt.US[jy].NAH2D.size(); hg++){
				Int_t thg = opt.US[jy].NAH2D[hg];
				for(Int_t ps = 0; ps<opt.PS[jz].NAH2D.size();ps++){
			    	    if(thg==opt.PS[jz].NAH2D[ps]){
			    		opt.Hyst2D[thg].HI.LieveTime+= 0.16*pam_events->GetTrigLevel2()->dltime[0];
					opt.Hyst2D[thg].HI.DeathTime+= 0.01*pam_events->GetTrigLevel2()->dltime[1];
				    }
				}
			    }
			    for(Int_t hg = 0; hg<opt.US[jy].NAH3D.size(); hg++){
				Int_t thg = opt.US[jy].NAH3D[hg];
				for(Int_t ps = 0; ps<opt.PS[jz].NAH3D.size();ps++){
			    	    if(thg==opt.PS[jz].NAH3D[ps]){
			    		opt.Hyst3D[thg].HI.LieveTime+= 0.16*pam_events->GetTrigLevel2()->dltime[0];
					opt.Hyst3D[thg].HI.DeathTime+= 0.01*pam_events->GetTrigLevel2()->dltime[1];
				    }
				}
			    }
			    jx = opt.US[jy].USP.size();
			}
		    }
	        }
	    }
	}
	//cout<<ntr<<endl;
	if(OneTrack){

	    if(ntr!=-1){
	    
	    	Bool_t verb = false;
	    	if(opt.opt.Nverb!=-1) if(i>=opt.opt.Nverb) verb = true;
		//if(opt.opt.Nverb==-1) verb = false;

	        if(verb) cout<<"\n\ni = "<<i<<"\n\n"<<endl;
	    			    	
	    	/*****************************************************************************/
    	    	// CALO Selection
    	    	/*****************************************************************************/
    	    	
    	    	Float_t caloqpre=track->GetCaloTrack()->qpre;
    	    	Int_t   calonpre=track->GetCaloTrack()->npre;
    	    	Float_t caloqtrack=track->GetCaloTrack()->qtrack;
    	    	Float_t caloqmax=pam_events->GetCaloLevel2()->qmax;
    	    	Int_t   dncore=track->GetCaloTrack()->ncore;
    	    	Int_t   lepton = 0;
    	    	Bool_t  ACSelection = false;
    	    	
    	    	if((hm0 & hs0)==0 && (hm1 & hs1)==0) ACSelection = true;
    	    	//cout<<"hm0 = "<<hm0<<" hm1 = "<<hm1<<" hs0 = "<<" hs1 = "<<hs1<<endl;
    	    	
    	    	Float_t betaev = 0.25*(track->GetToFTrack()->beta[0]+track->GetToFTrack()->beta[1]+track->GetToFTrack()->beta[2]+track->GetToFTrack()->beta[3]);
    	    	Float_t Ekin = 0;

    	    	if(detr<2. && Babs>0.24 && betaev>0.2){
    	    	    if(dncore > (pow(fabs(rigev),0.4)*700.-400.) && fabs(rigev)<11.2) lepton = 1;
    	    	    if(dncore > (pow(fabs(rigev),0.38)*700.-200.) && fabs(rigev)<48.) lepton = 1;
    	    	    if(dncore > 3200. && fabs(rigev)>=48.) lepton = 1;
    	    	    if(caloqtrack < (115.*fabs(rigev)-50.)) lepton = 0;
    	    	    if(dncore < (700.*pow(fabs(rigev),0.33)-400.) && fabs(rigev)<=70. && fabs(rigev)>32.) lepton = 0;
    	    	    if(fabs(rigev)>70. && dncore < 3200.) lepton = 0;
    	    	    if(dncore < (pow(fabs(rigev),0.38)*700.-200.) && fabs(rigev)>18. && fabs(rigev)<32.) lepton = 0;
    	    	    if(dncore < (700.*pow(fabs(rigev),0.4)-400.) && fabs(rigev)<18.) lepton = 0;
    	    	    if(fabs(rigev)<=40. && calonpre < pow(fabs(rigev),0.6)) lepton = 0;
    	    	    if(fabs(rigev)<=30. && calonpre > (70.-pow(fabs(rigev),0.4))) lepton = 0;
    	    	    if(fabs(rigev)>40. && calonpre < 10.) lepton = 0;
    	    	    if(fabs(rigev)>30. && calonpre > 32.) lepton = 0;
    	    	    if(caloqpre/caloqtrack > (exp(-fabs(rigev)/5.)+exp(-fabs(rigev)/0.1)+0.05)) lepton = 0;
    	    	    if(caloqmax/caloqtrack > (0.035+0.1/fabs(rigev))) lepton = 0;
    	    	}

		string Gname = "";
		//if(verb)for(Int_t o=0;o<opt.Hyst1D.size(); o++) cout<<"Name of Hyst1D["<<o<<"] is "<<opt.Hyst1D[o].HI.filename<<endl;
		//if(verb)for(Int_t o=0;o<opt.Hyst2D.size(); o++) cout<<"Name of Hyst2D["<<o<<"] is "<<opt.Hyst2D[o].HI.filename<<endl;
		
		if(verb){
		//if(rigev>1.7 && rigev<= 1.8 && L>6 && Lat>60 && Lat<-60){
		    cout<<"rigev = "<<rigev<<"\tdetr = "<<detr<<endl;
		    cout<<"Babs = "<<Babs<<"\tL = "<<L<<endl;
		}
		/*cout<<"Write Rigidity, please:"<<endl;
		cin>>rigev;
		cout<<"Write detr, please:"<<endl;
		cin>>detr;
		cout<<"Write lepton, please:"<<endl;
		cin>>lepton;
		cout<<"Write Anticounter, please:"<<endl;
		cin>>ACSelection;
		cout<<"Write Babs, please:"<<endl;
		cin>>Babs;
		cout<<"Write L, please:"<<endl;
		cin>>L;
		cout<<"Write Altitude, please:"<<endl;
		cin>>Alt;*/
		
		Double_t Rmin = 0;
		Double_t Rmax = 0;
		Double_t Detrmin = 0;
		Double_t Detrmax = 0;
		for(Int_t k=0;k<opt.SS.size();k++){
		    //Int_t NSS = opt.Hyst1D[k].NSS;
		    Int_t ntr = 0;
		    Int_t u = 0;
		    while(u<opt.SS[k].SP.size()){
		        Double_t Detrmin = 0;
		        Double_t Detrmax = 0;
		        //if(verb)cout<<"\t\trigev = "<<rigev<<endl;
		        //if(verb)cout<<"\t\tSelectionRmin["<<u<<"] = "<<opt.SS[k].SP[u].Rmin<<"\tSelectionRmax["<<u<<"] = "<<opt.SS[k].SP[u].Rmax<<endl;
		        if(rigev>opt.SS[k].SP[u].Rmin && rigev<=opt.SS[k].SP[u].Rmax){
		            //if(verb)cout<<"\t\tdetrminf["<<u<<"] = "<<opt.SS[k].SP[u].detrminf<<"\tdetrmaxf["<<u<<"] = "<<opt.SS[k].SP[u].detrmaxf<<endl;
		            if(opt.SS[k].SP[u].detrminf!=0){
			        switch ( opt.SS[k].SP[u].detrminf ) {
			    	    case 1:{Detrmin = 1/pow(fabs(rigev)+0.46,5)+1.1;if(verb){cout<<"\t\tCase 1: Detrmin = "<<Detrmin<<endl;}break;}
				    case 2:{Detrmin = 9.0*exp(-5.6*(fabs(rigev)-0.385))+3.7;if(verb){cout<<"\t\tCase 2: Detrmin = "<<Detrmin<<endl;}break;}
			    	    case 3:{Detrmin = 15*exp(-1.6*(fabs(rigev)-0.5))+4.;if(verb){cout<<"\t\tCase 3: Detrmin = "<<Detrmin<<endl;}break;}
				};
			    }else {Detrmin = opt.SS[k].SP[u].detrmore;if(verb){cout<<"\t\tDetrmore: Detrmin = "<<Detrmin<<endl;}}
			    if(opt.SS[k].SP[u].detrmaxf!=0){
			        switch ( opt.SS[k].SP[u].detrmaxf ) {
				    case 1:{Detrmax = 1/pow(fabs(rigev)+0.46,5)+1.1;if(verb){cout<<"\t\tCase 1: Detrmax = "<<Detrmax<<endl;}break;}
				    case 2:{Detrmax = 9.0*exp(-5.6*(fabs(rigev)-0.385))+3.7;if(verb){cout<<"\t\tCase 2: Detrmax = "<<Detrmax<<endl;}break;}
				    case 3:{Detrmax = 15*exp(-1.6*(fabs(rigev)-0.5))+4.;if(verb){cout<<"\t\tCase 3: Detrmax = "<<Detrmax<<endl;}break;}
				};
			    }else {Detrmax = opt.SS[k].SP[u].detrless;if(verb){cout<<"\t\tDetrless: Detrmax = "<<Detrmax<<endl;}}
			    Int_t Lept;
			    //if(verb)cout<<"\t\t Part Type is "<<opt.SS[k].PartType<<endl;
			    if(opt.SS[k].PartType==1 || opt.SS[k].PartType==4) Lept =1; else Lept = 0;
			    if(verb)cout<<"\t\tlepton = "<<lepton<<"\tLept = "<<Lept<<endl;
			    if(!(Detrmax==0 && Detrmin==0) && (detr<=Detrmin || detr>Detrmax)) ntr = -1;
			    //if(verb)cout<<"\t\tdetr = "<<detr<<"\tDetrmin = "<<Detrmin<<"\tDetrmax = "<<Detrmax<<endl;
			    if(opt.SS[k].SP[u].Calo && lepton==Lept) ntr=-1;
			    if(opt.SS[k].SP[u].AC && !ACSelection) ntr=-1;
			    //if(verb)cout<<"\t\tntr = "<<ntr<<endl;
			}
			u++;
		    }//while SP   
		    if(ntr!=-1) {
			vector<Int_t> N1D(0);
			vector<Int_t> N2D(0);
			Double_t Bmore = 0;
			Double_t Bless = 0;
			Double_t Lmore = 0;
			Double_t Lless = 0;
			Double_t altmore = 0;
			Double_t altless = 0;
			Double_t Rmore = 0;
			Double_t Rless = 0;
			Double_t Pmore = 0;
			Double_t Pless = 0;
			for(Int_t cu=0;cu<opt.SS[k].NAH1D.size();cu++){
			    Int_t nus=opt.Hyst1D[opt.SS[k].NAH1D[cu]].HI.NUS;
			    Int_t nps=opt.Hyst1D[opt.SS[k].NAH1D[cu]].HI.NPT;
			    for(Int_t jx = 0; jx<opt.US[nus].USP.size(); jx++){
				Bmore = opt.US[nus].USP[jx].Bmore;
			        Bless = opt.US[nus].USP[jx].Bless;
				Lmore = opt.US[nus].USP[jx].Lmore;
				Lless = opt.US[nus].USP[jx].Lless;
				altmore = opt.US[nus].USP[jx].altmore;
			        altless = opt.US[nus].USP[jx].altless;
				Rmore = opt.US[nus].USP[jx].Rmore;
			        Rless = opt.US[nus].USP[jx].Rless;
			        Pmore = opt.PS[nps].PitchMin;
			        Pless = opt.PS[nps].PitchMax;
				//cout<<"Pmore = "<<Pmore<<" Pless = "<<Pless<<endl;
				if(Babs>Bmore && Babs<=Bless && L>Lmore && L<=Lless && Alt>altmore && Alt<=altless && rigev>Rmore && rigev<=Rless && SB>Pmore && SB<=Pless){
			    	    //cout<<log10(rigev)<<endl;
			    	    N1D.resize(N1D.size()+1);
			    	    N1D[N1D.size()-1] = opt.SS[k].NAH1D[cu];
			    	    if(verb){
					cout<<"nus = "<<nus<<endl;
				        cout<<"Bmore = "<<Bmore<<"\tBless = "<<Bless<<endl;
					cout<<"Lmore = "<<Lmore<<"\tLless = "<<Lless<<endl;
				        cout<<"N1D["<<N1D.size()-1<<"] = "<<opt.SS[k].NAH1D[cu]<<endl;
					cout<<"HystName "<<opt.Hyst1D[N1D[N1D.size()-1]].HI.filename<<endl;
				        //Int_t we;
					//cin>>we;
				    }
			        jx = opt.US[nus].USP.size();
			        }
			    }
			}
			for(Int_t cu=0;cu<opt.SS[k].NAH2D.size();cu++){
			    Int_t nus=opt.Hyst2D[opt.SS[k].NAH2D[cu]].HI.NUS;
			    Int_t nps=opt.Hyst2D[opt.SS[k].NAH2D[cu]].HI.NPT;
			    for(Int_t jx = 0; jx<opt.US[nus].USP.size(); jx++){
				Bmore = opt.US[nus].USP[jx].Bmore;
			        Bless = opt.US[nus].USP[jx].Bless;
				Lmore = opt.US[nus].USP[jx].Lmore;
			        Lless = opt.US[nus].USP[jx].Lless;
				altmore = opt.US[nus].USP[jx].altmore;
			        altless = opt.US[nus].USP[jx].altless;
			        Rmore = opt.US[nus].USP[jx].Rmore;
			        Rless = opt.US[nus].USP[jx].Rless;
			        Pmore = opt.PS[nps].PitchMin;
			        Pless = opt.PS[nps].PitchMax;
				if(Babs>Bmore && Babs<Bless && L>Lmore && L<Lless && Alt>altmore && Alt<altless && rigev>Rmore && rigev<Rless && SB>Pmore && SB<Pless){
			    	    N2D.resize(N2D.size()+1);
			    	    N2D[N2D.size()-1] = opt.SS[k].NAH2D[cu];
			    	    if(verb){
			    	    //if(rigev>1.7 && rigev<= 1.8){
					cout<<"nus = "<<nus<<endl;
				        cout<<"Babs = "<<Babs<<" Bmore = "<<Bmore<<"\tBless = "<<Bless<<endl;
					cout<<"L = "<<L<<" Lmore = "<<Lmore<<"\tLless = "<<Lless<<endl;
				        cout<<"R = "<<rigev<<" Rmore = "<<Rmore<<"\tRless = "<<Rless<<endl;
				        cout<<"SB = "<<TMath::RadToDeg()*SB<<" Pmore = "<<Pmore<<"\tPless = "<<Pless<<endl;
				        cout<<"N2D["<<N2D.size()-1<<"] = "<<opt.SS[k].NAH2D[cu]<<endl;
				        cout<<"HystName "<<opt.Hyst2D[N2D[N2D.size()-1]].HI.filename<<endl;
					//Int_t we;
				        //cin>>we;
				    }
				jx = opt.US[nus].USP.size();
				}
			    }
			}
			for(Int_t d=0;d<N1D.size();d++){
			    Rmin = opt.ED[opt.Hyst1D[N1D[d]].HI.NED].Emin;
			    Rmax = opt.ED[opt.Hyst1D[N1D[d]].HI.NED].Emax;
			    Detrmin = opt.ED[opt.Hyst1D[N1D[d]].HI.NED].detrmin;
			    Detrmax = opt.ED[opt.Hyst1D[N1D[d]].HI.NED].detrmax;
			    Double_t Argv = 0;
			    if(rigev>Rmin && rigev<=Rmax && detr>Detrmin && detr<=Detrmax){
			        Int_t ExpType = opt.Hyst1D[N1D[d]].HI.ExpType;
			        Double_t LTI = 0;
			        Double_t DTI = 0;
			        if(ExpType == 1 || ExpType == 3) LTI = 0.16*pam_events->GetTrigLevel2()->dltime[0];
			        if(ExpType == 2 || ExpType == 3) DTI = 0.01*pam_events->GetTrigLevel2()->dltime[1];
			        switch  ( ExpType ){
			            case 1:{Argv = LTI;break;}
			            case 2:{Argv = DTI;break;}
				    case 3:{Argv = LTI+DTI;break;}
				    case 4:{
					    Argv = opt.Hyst1D[N1D[d]].HI.LieveTime;
					    opt.Hyst1D[N1D[d]].HI.LieveTime = 0;
					    break;
					    }
				    case 5:{
				            Argv = opt.Hyst1D[N1D[d]].HI.DeathTime;
					    opt.Hyst1D[N1D[d]].HI.DeathTime = 0;
					    break;
					    }
				    case 6:{
					    Argv = opt.Hyst1D[N1D[d]].HI.LieveTime+opt.Hyst1D[N1D[d]].HI.DeathTime;
					    opt.Hyst1D[N1D[d]].HI.LieveTime = 0;
					    opt.Hyst1D[N1D[d]].HI.DeathTime = 0;
					    break;
					    }
				};
				Double_t Xfill = 0;
				//cout<<"XType = "<<opt.Hyst1D[d].XType<<endl;
				//cout<<"Babs = "<<Babs<<endl;
				switch ( opt.Hyst1D[N1D[d]].XType ){
				    case 1:{ Xfill = TMath::RadToDeg()*PamAzim;break; }
				    case 2:{ Xfill = TMath::RadToDeg()*PamZenit;break; }
				    case 3:{ Xfill = SB; break; }
				    case 4:{ Xfill = detr; break; }
				    case 5:{ Xfill = rigev; break; }
				    case 6:{ Xfill = Lat; break; }
				    case 7:{ Xfill = Lon; break; }
				    case 8:{ Xfill = L; break; }
				    case 9:{ Xfill = Babs; break; }
				    case 10:{Xfill = Alt; break; }
				    case 11:{Xfill = LieveTime; break; }
				    case 12:{Xfill = DeathTime; break; }
				    case 13:{Xfill = SC; break; }
				    case 14:{Xfill = SA; break; }
				    case 17:{Xfill = log10(rigev); break;}
				};
				if(verb){
				    cout<<"rigev = "<<rigev<<"\tdetr = "<<detr<<endl;
				    cout<<"Babs = "<<Babs<<"\tL = "<<L<<endl;
				    cout<<"I Fill "<<opt.Hyst1D[N1D[d]].HI.filename<<" hystogramm"<<endl;
				    cout<<"Xfill = "<<Xfill<<" d = "<<d<<endl;
				}
				opt.Hyst1D[N1D[d]].Hyst1DF.Fill(Xfill);
				if(ExpType!=0) opt.Hyst1D[N1D[d]].Hyst1DFExp.Fill(Xfill,Argv);
				if(verb){
				    Int_t intt;
				    cin>>intt;
				}
			    }
			}//for d< N1D.size
			for(Int_t d=0;d<N2D.size();d++){
			    Rmin = opt.ED[opt.Hyst2D[N2D[d]].HI.NED].Emin;
			    Rmax = opt.ED[opt.Hyst2D[N2D[d]].HI.NED].Emax;
			    Detrmin = opt.ED[opt.Hyst2D[N2D[d]].HI.NED].detrmin;
			    Detrmax = opt.ED[opt.Hyst2D[N2D[d]].HI.NED].detrmax;
			    Double_t Argv = 0;
			    Double_t LTI = 0;
			    Double_t DTI = 0;
			    if(rigev>Rmin && rigev<=Rmax && detr>Detrmin && detr<=Detrmax){
			        Int_t ExpType = opt.Hyst2D[N2D[d]].HI.ExpType;
			        if(ExpType == 1 || ExpType == 3) LTI = 0.16*pam_events->GetTrigLevel2()->dltime[0];
			        if(ExpType == 2 || ExpType == 3) DTI = 0.01*pam_events->GetTrigLevel2()->dltime[1];
			        switch  ( ExpType ){
			            case 1:{Argv = LTI;break;}
			            case 2:{Argv = DTI;break;}
				    case 3:{Argv = LTI+DTI;break;}
				    case 4:{
					    Argv = opt.Hyst2D[N2D[d]].HI.LieveTime;
					    opt.Hyst2D[N2D[d]].HI.LieveTime = 0;
					    break;
					    }
				    case 5:{
				            Argv = opt.Hyst2D[N2D[d]].HI.DeathTime;
					    opt.Hyst2D[N2D[d]].HI.DeathTime = 0;
					    break;
					    }
				    case 6:{
					    Argv = opt.Hyst2D[N2D[d]].HI.LieveTime+opt.Hyst2D[N2D[d]].HI.DeathTime;
					    opt.Hyst2D[N2D[d]].HI.LieveTime = 0;
					    opt.Hyst2D[N2D[d]].HI.DeathTime = 0;
					    break;
					    }
				};
				Double_t Xfill = 0;
				Double_t Yfill = 0;
				//cout<<"YType = "<<opt.Hyst2D[d].YType<<endl;
				//cout<<"XType = "<<opt.Hyst2D[d].XType<<endl;
				switch ( opt.Hyst2D[N2D[d]].XType ){
				    case 1:{ Xfill = TMath::RadToDeg()*PamAzim;break; }
				    case 2:{ Xfill = TMath::RadToDeg()*PamZenit;break; }
				    case 3:{ Xfill = SB; break; }
				    case 4:{ Xfill = detr; break; }
				    case 5:{ Xfill = rigev; break; }
				    case 6:{ Xfill = Lat; break; }
				    case 7:{ Xfill = Lon; break; }
				    case 8:{ Xfill = L; break; }
				    case 9:{ Xfill = Babs; break; }
				    case 10:{Xfill = Alt; break; }
				    case 11:{Xfill = LieveTime; break; }
				    case 12:{Xfill = DeathTime; break; }
				    case 13:{Xfill = SC; break; }
				    case 14:{Xfill = SA; break; }
				    case 17:{Xfill = log10(rigev); break;}
				};
				switch ( opt.Hyst2D[N2D[d]].YType ){
				    case 1:{ Yfill = TMath::RadToDeg()*PamAzim;break; }
				    case 2:{ Yfill = TMath::RadToDeg()*PamZenit;break; }
				    case 3:{ Yfill = SB; break; }
				    case 4:{ Yfill = detr; break; }
				    case 5:{ Yfill = rigev; break; }
				    case 6:{ Yfill = Lat; break; }
				    case 7:{ Yfill = Lon; break; }
				    case 8:{ Yfill = L; break; }
				    case 9:{ Yfill = Babs; break; }
				    case 10:{Yfill = Alt; break; }
				    case 11:{Yfill = LieveTime; break; }
				    case 12:{Yfill = DeathTime; break; }
				    case 13:{Yfill = SC; break; }
				    case 14:{Yfill = SA; break; }
				    case 17:{Yfill = log10(rigev); break;}
				};
				if(verb){
				    cout<<endl;
				    cout<<"rigev = "<<rigev<<"\tdetr = "<<detr<<endl;
			    	    cout<<"Babs = "<<Babs<<"\tL = "<<L<<endl;
				    cout<<"I Fill "<<opt.Hyst2D[N2D[d]].HI.filename<<" hystogramm"<<endl;
				    cout<<"Xfill = "<<Xfill<<endl<<" Yfill = "<<Yfill<<endl;
				}
				opt.Hyst2D[N2D[d]].Hyst2DF.Fill(Xfill,Yfill);
				if(ExpType!=0) opt.Hyst2D[N2D[d]].Hyst2DFExp.Fill(Xfill,Yfill,Argv);
				if(verb){
				    Int_t intt;
				    cin>>intt;
				}
			    }
			}//for d<N"D.size
			
		    }//if ntr!=-1
		} // k<SS.size()
		if(verb)cout<<i<<endl;
	    } // if ntr!=-1;
	} // if GetNTrack==1;
    } //general loop
    
    TCanvas* Canv;
    TCanvas* Canvdiv;
    TCanvas* CanvExp;
    TCanvas* CanvLog;
    TString Ext[3];
    Ext[0] = ".jpg";	Ext[1] = ".png";	Ext[2] = ".bmp";
    
    if(opt.opt.verbose>=2) cout<<"\n\nSave block...\n\n";
    
    gStyle->SetPalette(1);
    
    for(Int_t i = 0; i<opt.Hyst1D.size(); i++){
        TH1F hdiv;
        TH1F h1n1D; TH1F h2n1D; TH1F hlog;
        if(opt.Hyst1D[i].HI.ExpType!=0){
    	    hlog.SetBins(80,-1.,3.);
    	    h1n1D.SetBins(opt.Hyst1D[i].Hyst1DF.GetNbinsX(),opt.Hyst1D[i].Hyst1DF.GetXaxis()->GetXmin(),opt.Hyst1D[i].Hyst1DF.GetXaxis()->GetXmax());
	    h2n1D.SetBins(opt.Hyst1D[i].Hyst1DF.GetNbinsX(),opt.Hyst1D[i].Hyst1DF.GetXaxis()->GetXmin(),opt.Hyst1D[i].Hyst1DF.GetXaxis()->GetXmax());
	    //h2div.SetBins(h12D->GetNbinsX(),h12D->GetXaxis()->GetXmin(),h12D->GetXaxis()->GetXmax(),h12D->GetNbinsY(),h12D->GetYaxis()->GetXmin(),h12D->GetYaxis()->GetXmax());
	    for(Int_t ii = 0; ii < opt.Hyst1D[i].Hyst1DF.GetNbinsX(); ii++){
		h1n1D.SetBinContent(ii,opt.Hyst1D[i].Hyst1DF.GetBinContent(ii)/opt.Hyst1D[i].Hyst1DF.GetEntries());
		h2n1D.SetBinContent(ii,opt.Hyst1D[i].Hyst1DFExp.GetBinContent(ii)/opt.Hyst1D[i].Hyst1DFExp.GetEntries());
	    }
    	    hdiv.SetName(opt.Hyst1D[i].Hyst1DF.GetName());
    	    hdiv.SetBins(opt.Hyst1D[i].Hyst1DF.GetNbinsX(),opt.Hyst1D[i].Hyst1DF.GetXaxis()->GetXmin(),opt.Hyst1D[i].Hyst1DF.GetXaxis()->GetXmax());
    	    hdiv.Divide(&h1n1D,&h2n1D);
    	    cout<<"i = "<<i<<endl;
    	    if(opt.Hyst1D[i].XType==17) hlog.Divide(&opt.Hyst1D[i].Hyst1DF,&binw);
        }
        if(opt.Hyst1D[i].HI.TypeExt > 0){
    	    Canv = new TCanvas("Canv",opt.Hyst1D[i].Hyst1DF.GetName(),1200,800);
    	    Canv->cd();
    	    opt.Hyst1D[i].Hyst1DF.Draw();
    	    Canv->SaveAs((TString)outdir+"ParDirCalc_"+opt.Hyst1D[i].HI.filename+".png");
    	    if(opt.Hyst1D[i].HI.ExpType!=0){
    		CanvExp = new TCanvas("CanvExp",opt.Hyst1D[i].Hyst1DFExp.GetName(),1200,800);
    		CanvExp->cd();
    		opt.Hyst1D[i].Hyst1DFExp.Draw();
    		CanvExp->SaveAs((TString)outdir+"ParDirCalc_"+opt.Hyst1D[i].HI.filename+"_Exp.png");
    		Canvdiv = new TCanvas("Canvdiv",hdiv.GetName(),1200,800);
    	        Canvdiv->cd();
    		hdiv.Draw();
    	        Canvdiv->SaveAs((TString)outdir+"ParDirCalc_"+opt.Hyst1D[i].HI.filename+"DivExp.png");
    	    }
    	    if(opt.Hyst1D[i].XType==17){
    		CanvLog = new TCanvas("CanvLog",opt.Hyst1D[i].Hyst1DF.GetName(),1200,800);
    		CanvLog->cd();
    	        hlog.Draw();
    		CanvLog->SaveAs((TString)outdir+"ParDirCalc_"+opt.Hyst1D[i].HI.filename+"_DivBinWide.png");
    	    }
    	}
    	if(opt.Hyst1D[i].HI.TypeExt == 0 || opt.Hyst1D[i].HI.TypeExt == 2){
    	    TFile fg((TString)outdir+"ParDirCalc_"+opt.Hyst1D[i].HI.filename+".root","recreate");
    	    opt.Hyst1D[i].Hyst1DF.Write();
    	    fg.Close();
    	    TFile fe((TString)outdir+"ParDirCalc_"+opt.Hyst1D[i].HI.filename+"_DivWide.root","recreate");
    	    hlog.Write();
    	    fe.Close();
    	    if(opt.Hyst1D[i].HI.ExpType!=0){
    		TFile fgExp((TString)outdir+"ParDirCalc_"+opt.Hyst1D[i].HI.filename+"_Exp.root","recreate");
    	        opt.Hyst1D[i].Hyst1DFExp.Write();
    		fg.Close();
    		TFile fgdiv((TString)outdir+"ParDirCalc_"+opt.Hyst1D[i].HI.filename+"_DivExp.root","recreate");
    	        hdiv.Write();
    		fgdiv.Close();
    	    }
    	    
        }
    }
    for(Int_t i = 0; i<opt.Hyst2D.size(); i++){
        TH2F hdiv;
        TH2F h1n2D; TH2F h2n2D;
        cout<<"i = "<<i<<endl;
        if(opt.Hyst2D[i].HI.ExpType!=0){
    	    //cout<<"Set hdiv"<<endl;
    	    h1n2D.SetBins(opt.Hyst2D[i].Hyst2DF.GetNbinsX(),opt.Hyst2D[i].Hyst2DF.GetXaxis()->GetXmin(),opt.Hyst2D[i].Hyst2DF.GetXaxis()->GetXmax(),opt.Hyst2D[i].Hyst2DF.GetNbinsY(),opt.Hyst2D[i].Hyst2DF.GetYaxis()->GetXmin(),opt.Hyst2D[i].Hyst2DF.GetYaxis()->GetXmax());
	    h2n2D.SetBins(opt.Hyst2D[i].Hyst2DF.GetNbinsX(),opt.Hyst2D[i].Hyst2DF.GetXaxis()->GetXmin(),opt.Hyst2D[i].Hyst2DF.GetXaxis()->GetXmax(),opt.Hyst2D[i].Hyst2DF.GetNbinsY(),opt.Hyst2D[i].Hyst2DF.GetYaxis()->GetXmin(),opt.Hyst2D[i].Hyst2DF.GetYaxis()->GetXmax());
	    //h2div.SetBins(h12D->GetNbinsX(),h12D->GetXaxis()->GetXmin(),h12D->GetXaxis()->GetXmax(),h12D->GetNbinsY(),h12D->GetYaxis()->GetXmin(),h12D->GetYaxis()->GetXmax());
	    for(Int_t ii = 0; ii < opt.Hyst2D[i].Hyst2DF.GetNbinsX(); ii++){
	        for(Int_t j = 0; j < opt.Hyst2D[i].Hyst2DF.GetNbinsY(); j++){
		    h1n2D.SetBinContent(ii,j,opt.Hyst2D[i].Hyst2DF.GetBinContent(ii,j)/opt.Hyst2D[i].Hyst2DF.GetEntries());
		    h2n2D.SetBinContent(ii,j,opt.Hyst2D[i].Hyst2DFExp.GetBinContent(ii,j)/opt.Hyst2D[i].Hyst2DFExp.GetEntries());
		}
	    }
    	    hdiv.SetName(opt.Hyst2D[i].Hyst2DF.GetName());
    	    hdiv.SetBins(opt.Hyst2D[i].Hyst2DF.GetNbinsX(),opt.Hyst2D[i].Hyst2DF.GetXaxis()->GetXmin(),opt.Hyst2D[i].Hyst2DF.GetXaxis()->GetXmax(),opt.Hyst2D[i].Hyst2DF.GetNbinsY(),opt.Hyst2D[i].Hyst2DF.GetYaxis()->GetXmin(),opt.Hyst2D[i].Hyst2DF.GetYaxis()->GetXmax());
    	    //cout<<"I try to divide it!!!"<<endl;
    	    hdiv.Divide(&h1n2D,&h2n2D);
    	    //cout<<"NbinsX1 = "<<h1n2D.GetNbinsX()<<" NBinsY1 = "<<h1n2D.GetNbinsY()<<endl;
    	    //cout<<"NbinsX2 = "<<h2n2D.GetNbinsX()<<" NBinsY2 = "<<h2n2D.GetNbinsY()<<endl;
    	    //cout<<"I has divide it!!! :))))\n0_o"<<endl;
        }
        //cout<<"opt.Hyst2D["<<i<<"].HI.TypeExt = "<<opt.Hyst2D[i].HI.TypeExt<<endl;
        if(opt.Hyst2D[i].HI.TypeExt > 0){
    	    //cout<<"opt.Hyst2D[i].HI.TypeExt = "<<opt.Hyst2D[i].HI.TypeExt<<endl;
    	    Canv = new TCanvas("Canv",opt.Hyst2D[i].Hyst2DF.GetName(),1200,800);
    	    Canv->cd();
    	    //cout<<"Canvas..."<<endl;
    	    opt.Hyst2D[i].Hyst2DF.Draw("colz");
    	    //cout<<"Draw hyst..."<<endl;
    	    Canv->SaveAs((TString)outdir+"ParDirCalc_"+opt.Hyst2D[i].HI.filename+".png");
    	    cout<<"Save canvas..."<<endl;
    	    if(opt.Hyst2D[i].HI.ExpType!=0){
    		CanvExp = new TCanvas("CanvExp",opt.Hyst2D[i].Hyst2DFExp.GetName(),1200,800);
    		CanvExp->cd();
    		opt.Hyst2D[i].Hyst2DFExp.Draw("colz");
    		CanvExp->SaveAs((TString)outdir+"ParDirCalc_"+opt.Hyst2D[i].HI.filename+"_Exp.png");
    	        Canvdiv = new TCanvas("Canvdiv",hdiv.GetName(),1200,800);
    		Canvdiv->cd();
    	        hdiv.Draw("colz");
    		Canvdiv->SaveAs((TString)outdir+"ParDirCalc_"+opt.Hyst2D[i].HI.filename+"DivExp.png");
    	    }
        }
        if(opt.Hyst2D[i].HI.TypeExt == 0 || opt.Hyst2D[i].HI.TypeExt == 2){
    	    TFile fg((TString)outdir+"ParDirCalc_"+opt.Hyst2D[i].HI.filename+".root","recreate");
    	    opt.Hyst2D[i].Hyst2DF.Write();
    	    fg.Close();
    	    if(opt.Hyst2D[i].HI.ExpType!=0){
    		TFile fgExp((TString)outdir+"ParDirCalc_"+opt.Hyst2D[i].HI.filename+"_Exp.root","recreate");
    	        opt.Hyst2D[i].Hyst2DFExp.Write();
    		fg.Close();
    		cout<<"I'm here!"<<endl;
    		TFile fgdiv((TString)outdir+"ParDirCalc_"+opt.Hyst2D[i].HI.filename+"_DivExp.root","recreate");
    	        hdiv.Write();
    		fgdiv.Close();
    		cout<<"I'm here!"<<endl;
    	    }
        }
        //delete Canv;
        //delete Canvdiv;
        //delete CanvExp;
        //delete CanvLog;
    }

    /**************************************************************************************/
    //Draw hystogramms for protons. Volodia
    /**************************************************************************************/

/*    
    plogCanvas->cd();
    //plogCanvas->SetLogx();
    plogCanvas->SetLogy();
    plogCanvas->SetGrid(); 
    //BinLogX(proton1);
    Float_t Energybin[80],Fluxbin[80],dE[80],dFlux[80];

    for (Int_t i=0;i<80;i++){
	dFlux[i]=proton1log->GetBinContent(i+1);
	dFlux[i]=1./sqrt(dFlux[i]);
	// bintime3eq[i]=time3eq->GetBinContent(i+1);	   
    }

    proton1log->Divide(proton1log,binw);
    proton1log->Divide(proton1log,Geomfactorlog);
    proton1log->Scale(0.001);  //MeV ->GeV
    proton1log->Scale(10000.); //cm2 -> m2
    //proton1log->Scale(1./exposure[0]);

    for (Int_t i=0;i<80;i++){
	Fluxbin[i]=proton1log->GetBinContent(i+1);
        dFlux[i]=Fluxbin[i]*dFlux[i];
	Energybin[i]=pow(10.,(Float_t)i*0.05+0.025-1);
        Energybin[i+1]=pow(10.,(Float_t)(i+1)*0.05+0.025-1.);
	dE[i]=(Energybin[i+1]-Energybin[i])/2.;
	// bintime3eq[i]=time3eq->GetBinContent(i+1);
	//cout<<i<<"  "<<bintime3[i]<<"  "<<bintime3eq[i]<<endl;
	//flux_out<<Energybin[i]<<"  "<<Fluxbin[i]<<" "<<dE[i]<<" "<<dFlux[i]<<endl;
    }

    proton1log->Draw("");

*/
}