--- ParDirCalc/QtoInclination.C	2008/01/17 16:15:13	1.1.1.1
+++ ParDirCalc/QtoInclination.C	2008/09/16 09:38:49	1.3
@@ -6,6 +6,7 @@
 //
 #include "TrAng.h"
 #include "QtoInc.h"
+#include "option.h"
 
 #include <fstream>
 #include <iostream>
@@ -24,15 +25,370 @@
 #include <TCanvas.h>
 //#include <TMatrixD.h>
 #include "TStyle.h"
-#include <TH2F.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;
@@ -41,13 +397,39 @@
 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")){
-    cout<<"\nUsage \n";
+    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";
@@ -61,33 +443,590 @@
     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[3],"--OutPuth")){
-    if (4 >= argc){
-    cout<<"--OutPath needs argument\n See --help\n";    
-    }else outdir = argv[4];
+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);
 }
 
-for(Int_t i=4; i<argc; i++){
-    if(!strcmp(argv[i],"-DoTr")) DoTr = true;
-    if(!strcmp(argv[i],"-PamEff")) PamEff = true;
-    if(!strcmp(argv[i],"-DiffHist")){DiffHist = true; DoTr = true;}
-    if(!strcmp(argv[i],"-AngHist")) AnglHist = true;
-    if(!strcmp(argv[i],"-PamAngTime")) PamAngTime = true;
-    if(!strcmp(argv[i],"-PamExp")) PamExp = true;
+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(!PamEff && !DiffHist && !AnglHist && !PamAngTime && !PamExp) AnglHist = true;
+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";
+	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);
-    }else Calculate(argv[2],outdir,AnglHist,DiffHist,PamEff,DoTr,PamAngTime,PamExp);
 }
 
+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
@@ -116,173 +1055,90 @@
  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){
-    
+
+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
     /**********************************************/
     
-    //gStyle->SetOptStat(111111);
-    TSystemDirectory *workdir = new TSystemDirectory("workdir",dirname);
-    TList *flist=workdir->GetListOfFiles();
+    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(); 
+    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";
+    
     
-    cout<<"Number of events: "<<nevents<<endl;
     /********************************************************************************/
     /*****NEED TO  CHANGE FOR OTHER COMPUTERS****************************************/
-    pam_events->GetTrkLevel2()->LoadField("/data01/pamhome/installed/pamela_software/calib/trk-param/field_param-0/");
+    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;
-    Double_t Py;
-    Double_t Pz;
-    
-    /*********************************************************************************/
-    // Angular exposure. How much time Pamela have any Azimutal and Zenital angle
-    /*********************************************************************************/
-    
-    TH2F *hhigh;
-    TH2F *hlow;
-    
-    TFile fhigh((TString)outdir+"PamAngEfficiencyHighEnergy.root");
-    TFile flow((TString)outdir+"PamAngEfficiencyLowEnergy.root");
-    if(PamExp){
-    
-	if (fhigh.IsZombie()||flow.IsZombie()){
-	    cout<<"Problem with Hystogrammfiles"<<endl;
-	    exit(1);
-	}else{
-	    hhigh = (TH2F*)fhigh.Get("PamAngEffhigh");
-	    hlow = (TH2F*)flow.Get("PamAngEfflow");
-	}
-    
-    }
+    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(AnglHist){
-    
-    TH1F *PitchExpositionWithoutBrazil = new TH1F("PitchExpositionWithoutBrazil", "Pitch Exposition without Brazil with Track", 360, 0, 180);
-    TH1F *PitchExpositionEquator = new TH1F("PitchExpositionEquator", "Pitch Exposition in Equator  with Track", 360, 0, 180);
-    TH1F *PitchExpositionBrazil = new TH1F("PitchExpositionBrazil", "Pitch Exposition Brazil with Track", 360, 0, 180);
-    TH1F *MainAxisPamelaExpositionWithoutBrazil = new TH1F("MainAxisPamelaPitchExpositionWithoutBrazil", "Main Axis of Pamela Exposition without Brazil with Track", 360, 0, 180);
-    TH1F *MainAxisPamelaExpositionEquator = new TH1F("MainAxisPamelaPitchExpositionEquator", "Main Axis of Pamela Exposition in Equator with Track", 360, 0, 180);
-    TH1F *MainAxisPamelaExpositionBrazil = new TH1F("MainAxisPamelaExpositionBrazil", "Main Axis of Pamela Exposition Brazil with Track", 360, 0, 180);
-    TH1F *PitchExpositionNoOrientationWithoutBrazil = new TH1F("PitchExpositionNoOrientationWithoutBrazil", "Pitch Exposition without orientation Information without Brazil with Track", 360, 0, 180);
-    TH1F *PitchExpositionNoOrientationEquator = new TH1F("PitchExpositionNoOrientationEquator", "Pitch Exposition without orientation Information in Equator with Track", 360, 0, 180);
-    TH1F *PitchExpositionNoOrientationBrazil = new TH1F("PitchExpositionNoOrientationBrazil", "Pitch Exposition without orientation Information in Brazil with Track", 360, 0, 180);
-    TH1F *DiferenceBetweenNoOrientationAndWithOrientationWithoutBrazil = new TH1F("DiferenceBetweenNoOrientationAndWithOrientationWithoutBrazil", "Diference between no orientation data and with orientation without Brazil with Track", 360, 0, 180);
-    TH1F *DiferenceBetweenNoOrientationAndWithOrientationEquator = new TH1F("DiferenceBetweenNoOrientationAndWithOrientationEquator", "Diference between no orientation data and with orientation in Equator with Track", 360, 0, 180);
-    TH1F *DiferenceBetweenNoOrientationAndWithOrientationBrazil = new TH1F("DiferenceBetweenNoOrientationAndWithOrientationBrazil", "Diference between no orientation data and with orientation in Brazil with Track", 360, 0, 180);
-    TH1F *DifferenceCountWithoutBrazil = new TH1F("DiferenceCountWithoutBrazil","Diference without Brazil",150,0.001,3);
-    TH1F *DifferenceCountEquator = new TH1F("DiferenceCountWithoutBrazil","Diference without Brazil",150,0.001,3);
-    TH1F *DifferenceCountInBrazil = new TH1F("DiferenceCountInBrazil","Diference in Brazil",150,0.001,3);
-    
-    TH1F *PitchCountWithoutBrazil = new TH1F("PitchCountWithoutBrazil", "Pitch Count without Brazil with Track", 360, 0, 180);
-    TH1F *PitchCountEquator = new TH1F("PitchCountEquator", "Pitch Count in Equator with Track", 360, 0, 180);
-    TH1F *PitchCountBrazil = new TH1F("PitchCountBrazil", "Pitch Count Brazil with Track", 360, 0, 180);
-    TH1F *MainAxisPamelaCountWithoutBrazil = new TH1F("MainAxisPamelaPitchCountWithoutBrazil", "Main Axis of Pamela Count without Brazil with Track", 360, 0, 180);
-    TH1F *MainAxisPamelaCountEquator = new TH1F("MainAxisPamelaPitchCountEquator", "Main Axis of Pamela Count in Equator with Track", 360, 0, 180);
-    TH1F *MainAxisPamelaCountBrazil = new TH1F("MainAxisPamelaCountBrazil", "Main Axis of Pamela Count Brazil with Track", 360, 0, 180);
-    TH1F *PitchCountNoOrientationWithoutBrazil = new TH1F("PitchCountNoOrientationWithoutBrazil", "Pitch Count without orientation Information without Brazil with Track", 360, 0, 180);
-    TH1F *PitchCountNoOrientationEquator = new TH1F("PitchCountNoOrientationEquator", "Pitch Count without orientation Information in Equator with Track", 360, 0, 180);
-    TH1F *PitchCountNoOrientationBrazil = new TH1F("PitchCountNoOrientationBrazil", "Pitch Count without orientation Information in Brazil with Track", 360, 0, 180);
-    //TH1F *Diference2BetweenNoOrientationAndWithOrientationWithoutBrazil = new TH1F("Diference2BetweenNoOrientationAndWithOrientationWithoutBrazil", "Diference between no orientation data and with orientation without Brazil with Track", 360, 0, 180);
-    //TH1F *Diference2BetweenNoOrientationAndWithOrientationEquator = new TH1F("Diference2BetweenNoOrientationAndWithOrientationEquator", "Diference between no orientation data and with orientation in Equator with Track", 360, 0, 180);
-    //TH1F *Diference2BetweenNoOrientationAndWithOrientationBrazil = new TH1F("Diference2BetweenNoOrientationAndWithOrientationBrazil", "Diference between no orientation data and with orientation in Brazil with Track", 360, 0, 180);
-    //TH1F *Difference2CountWithoutBrazil = new TH1F("DiferenceCountWithoutBrazil","Diference without Brazil",150,0.001,3);
-    //TH1F *Difference2CountEquator = new TH1F("DiferenceCountWithoutBrazil","Diference without Brazil",150,0.001,3);
-    //TH1F *Difference2CountInBrazil = new TH1F("DiferenceCountInBrazil","Diference in Brazil",150,0.001,3);
-    //TH1F *EastPitchAllCount = new TH1F("EastPitchAllCount","East Count",)
-    
-    TH1F *ZenitAngle = new TH1F("ZenitAngle","Zenit Angle",360,120,180);
-    //TH1F *ZenitAngleEquator = new TH1F("ZenitAngleWithoutEquator","Zenit Angle in Equator",360,120,180);
-    //TH1F *ZenitAngleBrazil = new TH1F("ZenitAngleBrazil","Zenit Angle in Brazil",360,120,180);
-    
-    TH2F *PitchvsLatInBrazil = new TH2F("PitchvsLatInBrazil","Pitch Angle, deg; Latitude, deg",360,-40,-5,360,40,110);
-    TH2F *LieveTimevsLatInBrazil = new TH2F("LieveTimevsLatInBrazil","Live time, deg; Latitude deg",360,-40,-5,1000,0,80);
-    
-    
-    TCanvas *c1 = new TCanvas("c1","Exposition and count in Brazil",1200,3000);
-    TCanvas *c2 = new TCanvas("c2","Exposition and Count without Brazil",1200,3000);
-    TCanvas *c3 = new TCanvas("c3","Exposition and Count on Equator",1200,3000);
-    TCanvas *c4 = new TCanvas("c4","Pitch-Angle exposition and count",1200,3400);
-    TCanvas *c5 = new TCanvas("c5","PitchAndLieveTimevsLatitude",1200,800);
-    TCanvas *c6 = new TCanvas("c6","Diferense",1200,800);
-    //TCanvas *c6 = new TCanvas("c6","CountEquator",1200,3000);
-    
-    //}
+    if(opt.opt.verbose>=2) cout<<"\n--AnglHyst Hystogramm block...\n\n";
     
     /***************************************************************************************/
-    //Histogramms for diferences in results between calculating using DoTrack2
+    //Histogramms for differences in results between calculating using DoTrack2
     //and using axv[0], ayv[0]
     /**************************************************************************************/
     
-    //if(DiffHist){
-    
-    TH1F *DifDoTr2ANdaxvAzim = new TH1F("DifDoTr2AndaxvAzim","Diference between Azimutal angles",360,-180,180);
-    TH1F *DifDoTr2ANdaxvZenit = new TH1F("DifDoTr2AndaxvZenit","Diference between Zenit angles",360,-180,180);
-    
-    //}
-    
-    /*************************************************************************************/
-    //Histogramm for Pamela's angular efficiency
-    /*************************************************************************************/
-    
-    TH2F *PamAngEffhigh = new TH2F("PamAngEffhigh","Zenit Angle, deg; Azimutal Angle, deg",360,0,360,100,80,180);
-    TH2F *PamAngEfflow = new TH2F("PamAngEfflow","Zenit Angle, deg; Azimutal Angle, deg",360,0,360,100,80,180);
-    TCanvas *c7 = new TCanvas("c7","Pamela's angular efficiency", 1600,2000);
-    
-    /*************************************************************************************/
-    //Hystogramm to count how much time Pamela have each PitchAngle
-    /*************************************************************************************/
-    
-    TH1F *PitchTime = new TH1F("PitchTime","Time for Pitch",360,0,360);
-    TCanvas *c8 = new TCanvas("c8","PitchTime",1200,600);
-    
-    TH1F *PitchExposure = new TH1F("PitchExposure","Time for Pitch",360,0,360);
-    TCanvas *c9 = new TCanvas("c9","PitchTime",1200,600);
+    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.);
-    proton1log->GetXaxis()->SetTitle("log10(Energy) , GeV");
-    proton1log->SetLineColor(kRed);
-    TCanvas *plogCanvas = new TCanvas("proton flux","protonflux", 800, 800);
+    
+//    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];
+    //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.)); 
+	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);
-    }
-*/    
+    //for (Int_t l=0 ; l<80; l++) {
+    //binwide[l]=binw->GetBinContent(l+1);
+//    }
+
     /****************************************************************************/
     //Geometrical Factor. Volodia
     /****************************************************************************/
@@ -305,458 +1161,662 @@
     //General loop
     /****************************************************************************/
 
-    for(ULong_t i=0; i<nevents; i++){
+    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;
 	
-	if(pam_events->GetTrkLevel2()->GetNTracks()==1){
-	    
-	    /*************************************************************************/
-	    //Getting general parameters
-	    /*************************************************************************/
-	    
-	    track = pam_events->GetTrack(0);
-	    Int_t M0DE = pam_events->GetOrbitalInfo()->mode;  //0 is zero here
-	    Float_t Bx = -pam_events->GetOrbitalInfo()->Bdown;
-	    Float_t By = pam_events->GetOrbitalInfo()->Beast;
-	    Float_t Bz = pam_events->GetOrbitalInfo()->Bnorth;
-	    Float_t Babs = pam_events->GetOrbitalInfo()->Babs;
-	    Float_t L = pam_events->GetOrbitalInfo()->L;
-	    Float_t rigev = 1/track->GetTrkTrack()->GetDeflection();
-	    
-	    /*************************************************************************/
-	    //Track selection
-	    /*************************************************************************/
-	    
-	    ntr = 0;
-	    if ((track->GetTrkTrack()->chi2<=0) || (track->GetTrkTrack()->chi2>=500)) ntr=-1;
-	    for(Int_t ii=1;ii<=12;ii++){
-		if ((track->GetToFTrack()->beta[i]<0) || (track->GetToFTrack()->beta[i]>99)) ntr=-1;
-	    }
-	    if((track->GetTrkTrack()->GetNX()<=4)&&(track->GetTrkTrack()->GetNY()<4)) ntr=-1;
-	    if((M0DE!=0)&&(M0DE!=1)&&(M0DE!=6)&&(M0DE!=2)&&(M0DE!=3)&&(M0DE!=8)&&(M0DE!=4)) ntr=-1;
-	    if(sqrt(pow(Bx,2)+pow(By,2)+pow(Bz,2))<0.001) ntr=-1;
-	    if(rigev<0.600) ntr=-1;
-	    
-	    if(ntr!=-1 || PamExp){
-	    
-		/*******************************************************************************/
-		//Transit reference system
-	        /*******************************************************************************/
-	    
-		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);
-	    
-	        /*****************************************************************************/
-		//Calculate Zenit and Azimutal angle and vector of particle flight in Pamela
-	        //using axv[0] and ayv[0] variables
-		/*****************************************************************************/	    
+	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(!DoTr || DiffHist){
+	    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();
+	        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);
-	    
+	    	
+	    	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
-		/*****************************************************************************/
+	    /*****************************************************************************/
+	    //Calculate Zenit and Azimutal angle and vector of particle flight in Pamela
+	    //using DoTrack2 procedure
+	    /*****************************************************************************/
 	    
-		if(DoTr){
+	    if(opt.opt.DoTr){
 	    
-	        track->GetTrkTrack()->DoTrack2(tr);
-		Double_t E11x = tr->x[0];
+		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 E22x = tr->x[3];
 	        Double_t E22y = tr->y[3];
-		Double_t E22z = zin[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 = TMath::RadToDeg()*atan(TMath::Abs(E22y-E11y)/TMath::Abs(E22x-E11x));
 	        MyAzim = MyAzim;//-180;
-		MyZenit = TMath::RadToDeg()*acos((E22z-E11z)/norm);
+	        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 = MyAzim;
 	        if(E22x-E11x<0 && E22y-E11y >0) MyAzim = 180. - 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;
+	        Py = (E22y-E11y)/norm;
 	        Pz = (E22z-E11z)/norm;
-	    
-		}
-	    
-		/***************************************************************************************/
-	        //Fill histogramms for diferences in results between calculating using DoTrack2
-		//and using axv[0], ayv[0]
-	        /**************************************************************************************/
-	    
-		if(DiffHist){
-	    
-		DifDoTr2ANdaxvZenit->Fill(MyZenit-(TMath::RadToDeg()*PamZenit));
-		DifDoTr2ANdaxvAzim->Fill(MyAzim-(TMath::RadToDeg()*PamAzim));
-	    
-	        }
-	    	    
-	        if(PamEff||AnglHist||PamExp){
-
-		/***********************************************************************************/
-	        //Transit vector in Pamela reference system to GEO
-		/***********************************************************************************/
-	    
-		TMatrixD Eij = PO->PamelatoGEO(Iij,Px,Py,Pz);
-	        Double_t A1 = Iij(0,2);
-		Double_t A2 = Iij(1,2);
-	        Double_t A3 = Iij(2,2);
-
-		/*******************************************************************************/
-		//Calculating angles (Pitch, Pamela's main axis etc.)
-		/*******************************************************************************/
-		
-		Double_t SB = PO->GetPitchAngle(Eij(0,0),Eij(1,0),Eij(2,0),Bx,By,Bz); 	// Pitch angle
-		Double_t SA = PO->GetPitchAngle(A1,A2,A3,Bx,By,Bz);			// Angle between Pamela's main axiz and B
-	        Double_t SC = PO->GetPitchAngle(1,0,0,Bx,By,Bz);			// Angle between direction to xenit and B
-	        Double_t ZC = PO->GetPitchAngle(Px,Py,Pz,0,0,1);			// Zenit Angle of Particle flight
 	        
-		/******************************************************************************/
-		//Proton selection. Volodia.
-		/******************************************************************************/
-/*		
-		if(14.9/L/L > 0. && 14.9/L/L < 1. & Babs >0.24) {
-    		Int_t detr = 0;
-    		//if((rigev < -0.8  && betaev >0.8)||(rigev>-0.8 && rigev < 0.&& betaev > 0.8) ) //electron2->Fill(log10(fabs(rigev)));
-		    if(detr  < (15*exp(-1.6*(fabs(rigev)-0.5))+4)  && rigev >0.5  || (detr < 30 && fabs(rigev) <0.5)) {
-			Float_t Ekin=ConvertR2T(rigev,0.938,1., 1.);
-			proton1log->Fill(log10(Ekin));
-		    }
-		}
-*/		
-		/************/
-		//Exposure
-		/************/
-		
-		//for (Int_t itime=0; itime<16;itime++){
-		//    if(14.9/L/L > itime && 14.9/L/L < (itime+1) & Babs >0.24) exposure[itime]=exposure[itime]+livetime/1000.;}
-		
-		/***********************************************************************/
-		//Polar area. To calculate angular efficiencÐÐy of registration
-		/***********************************************************************/
-		
-		if(PamExp){
-
-		Float_t AnglTime = 0.16*pam_events->GetTrigLevel2()->dltime[0]+0.01*pam_events->GetTrigLevel2()->dltime[1];
-		for(Int_t i=0;i<=360;i++){
-		    for(Int_t j=0;j<=100;j++){
-		        if(hhigh->GetBin(i,j)>0){
-		    	    Double_t pamZenit=(j+80-0.5)*TMath::DegToRad();
-		    	    Double_t pamAzim=(i*-0.5)*TMath::DegToRad();
-		    	    Double_t _z = cos(pamZenit);
-			    Double_t _x = sin(pamZenit)*cos(pamAzim);
-		    	    Double_t _y = sin(pamZenit)*sin(pamAzim);
-			    TMatrixD Uij = PO->PamelatoGEO(Iij,_x,_y,_z);
-		    	    Double_t SP = PO->GetPitchAngle(Uij(0,0),Uij(1,0),Uij(2,0),Bx,By,Bz);
-			    PitchExposure->Fill(SP,AnglTime);
- 			}
-		    }
-		}
-		
-		}
-		
-		if(PamEff){
+	    }
+	    	
+	    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);
 		
-		if(L>6){
-		Float_t Mass = 0;
-    		Int_t detr = 0;
-    		Float_t Ekin = 0;
-    		Float_t betaev = 0.25*(track->GetToFTrack()->beta[0]+track->GetToFTrack()->beta[1]+track->GetToFTrack()->beta[2]+track->GetToFTrack()->beta[3]);    
-		    if(14.9/L/L > 0. && 14.9/L/L < 1. & Babs >0.24) {	
-    			if((rigev < -0.8  && betaev >0.8)||(rigev>-0.8 && rigev < 0.&& betaev > 0.8) ) {
-    			    Ekin=ConvertR2T(rigev,0.511,-1.);
-    			    if (Ekin>1) PamAngEffhigh->Fill(TMath::RadToDeg()*PamAzim,TMath::RadToDeg()*PamZenit);
-    			    if (Ekin<1) PamAngEfflow->Fill(TMath::RadToDeg()*PamAzim,TMath::RadToDeg()*PamZenit);
-    			}//electron2->Fill(log10(fabs(rigev)));
-			
-			if(detr  < (15*exp(-1.6*(fabs(rigev)-0.5))+4)  && rigev >0.5  || (detr < 30 && fabs(rigev) <0.5)) {
-			    Ekin=ConvertR2T(rigev,0.938,-1.);
-			    if (Ekin>1) PamAngEffhigh->Fill(TMath::RadToDeg()*PamAzim,TMath::RadToDeg()*PamZenit);
-			    if (Ekin<1) PamAngEfflow->Fill(TMath::RadToDeg()*PamAzim,TMath::RadToDeg()*PamZenit);
+	    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();
 			}
 		    }
-		}
-		
-		}
-		
-		/***********************************************************************/
-		//Filling hystogram to count how much time Pamela have each Pitch-Angle
-		/***********************************************************************/
-		
-		if(PamAngTime) {
-
-		Float_t AnglTime = 0.16*pam_events->GetTrigLevel2()->dltime[0]+0.01*pam_events->GetTrigLevel2()->dltime[1];
-		PitchTime->Fill(SA,AnglTime);
-
-		}
-
-		if(AnglHist){
+	        }
+	    }
+	}
+	//cout<<ntr<<endl;
+	if(OneTrack){
 
-		/***********************************************************************/
-		//Filling angular histogramms for Brazil
-		/***********************************************************************/
-
-		if((Babs<0.19)&&(L<1.75)){
-			Argv = 0.16*pam_events->GetTrigLevel2()->dltime[0];
-			PitchExpositionBrazil->Fill(SB,Argv);
-			MainAxisPamelaExpositionBrazil->Fill(SA,Argv);
-			PitchExpositionNoOrientationBrazil->Fill(SC,Argv);
-			PitchCountBrazil->Fill(SB);
-			MainAxisPamelaCountBrazil->Fill(SA);
-			PitchCountNoOrientationBrazil->Fill(SC);
-			ZenitAngle->Fill(ZC);
-			PitchvsLatInBrazil->Fill(pam_events->GetOrbitalInfo()->lat,SA);
-			LieveTimevsLatInBrazil->Fill(pam_events->GetOrbitalInfo()->lat,Argv);
-		}
-		
-		/***********************************************************************/
-		//Filling angulars histogramms for areas without Brazil
-		/***********************************************************************/
+	    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(Babs>0.24){
-			Argv = 0.16*pam_events->GetTrigLevel2()->dltime[0];
-			PitchExpositionWithoutBrazil->Fill(SB,Argv);
-			MainAxisPamelaExpositionWithoutBrazil->Fill(SA,Argv);
-			PitchExpositionNoOrientationWithoutBrazil->Fill(SC,Argv);
-			PitchCountWithoutBrazil->Fill(SB);
-			MainAxisPamelaCountWithoutBrazil->Fill(SA);
-			PitchCountNoOrientationWithoutBrazil->Fill(SC);
-			ZenitAngle->Fill(ZC);
-			
-			/***************************************************************/
-			//Filling angulars histogramm for equator
-			/***************************************************************/
-			
-			if(L<1.2){
-			    Argv = 0.16*pam_events->GetTrigLevel2()->dltime[0];
-			    PitchExpositionEquator->Fill(SB,Argv);
-			    MainAxisPamelaExpositionEquator->Fill(SA,Argv);
-			    PitchExpositionNoOrientationEquator->Fill(SC,Argv);
-			    PitchCountEquator->Fill(SB);
-			    MainAxisPamelaCountEquator->Fill(SA);
-			    PitchCountNoOrientationEquator->Fill(SC);
-			    ZenitAngle->Fill(ZC);
-			}
+		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;*/
 		
-		} // if(AnglHist)
-	    
-	    } // if(PamEff||AnglHist)
-	    
+		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
-
-    Double_t BinContent;
     
-    /*********************************************************/
-    //Divide canvases for angular histogramms
-    /*********************************************************/
-    
-    if(AnglHist){
-    
-    c1->Divide(1,6);
-    c2->Divide(1,6);
-    c3->Divide(1,6);
-    c4->Divide(1,7);
-    c5->Divide(1,2);
+    TCanvas* Canv;
+    TCanvas* Canvdiv;
+    TCanvas* CanvExp;
+    TCanvas* CanvLog;
+    TString Ext[3];
+    Ext[0] = ".jpg";	Ext[1] = ".png";	Ext[2] = ".bmp";
     
-    }
-    
-    /***************************************************************************************/
-    //Canvas for histogramms for diferences in results between calculating using DoTrack2
-    //and using axv[0], ayv[0]
-    /**************************************************************************************/
-    
-    if(DiffHist){
-    
-    c6->Divide(1,2);    
-    
-    }
-    
-    /**************************************************************************************/
-    //Canvas for Pamela's efficiency hystogramm
-    /**************************************************************************************/
+    if(opt.opt.verbose>=2) cout<<"\n\nSave block...\n\n";
     
-    if(PamEff){
-    
-    c7->Divide(1,2);
-    
-    }
+    gStyle->SetPalette(1);
     
-    /**************************************************************************************/
-    //Draw Angular Histogramms
-    /**************************************************************************************/
-
-    if(AnglHist){
-
-    c1->cd(1);
-    MainAxisPamelaExpositionBrazil->Draw();
-    c1->cd(2);
-    MainAxisPamelaCountBrazil->Draw();
-    c1->cd(3);
-    PitchExpositionNoOrientationBrazil->Draw();
-    c1->cd(4);
-    PitchCountNoOrientationBrazil->Draw();
-    c1->cd(5);
-    DiferenceBetweenNoOrientationAndWithOrientationBrazil->Divide(PitchExpositionNoOrientationBrazil,MainAxisPamelaExpositionBrazil);
-    DiferenceBetweenNoOrientationAndWithOrientationBrazil->Draw();
-    c1->cd(6);
-    for(Int_t iu=0;iu<360;iu++){
-	BinContent = DiferenceBetweenNoOrientationAndWithOrientationBrazil->GetBinContent(iu);
-	DifferenceCountInBrazil->Fill(BinContent);
-    }
-    DifferenceCountInBrazil->Draw();
-    //DifferenceCountInBrazil->SaveAs((TString)outdir+"PIC001A_BrazilExpositionAndCount.root");
-    c1->SaveAs((TString)outdir+"PIC001A_BrazilExpositionAndCount.jpg");
-    
-    c2->cd(1);
-    MainAxisPamelaExpositionWithoutBrazil->Draw();
-    c2->cd(2);
-    MainAxisPamelaCountWithoutBrazil->Draw();
-    c2->cd(3);
-    PitchExpositionNoOrientationWithoutBrazil->Draw();
-    c2->cd(4);
-    PitchCountNoOrientationWithoutBrazil->Draw();
-    c2->cd(5);
-    DiferenceBetweenNoOrientationAndWithOrientationWithoutBrazil->Divide(PitchExpositionNoOrientationWithoutBrazil,MainAxisPamelaExpositionWithoutBrazil);
-    DiferenceBetweenNoOrientationAndWithOrientationWithoutBrazil->Draw();
-    c2->cd(6);
-    for(Int_t iu=0;iu<360;iu++){
-	BinContent = DiferenceBetweenNoOrientationAndWithOrientationWithoutBrazil->GetBinContent(iu);
-	DifferenceCountWithoutBrazil->Fill(BinContent);
-    }
-    DifferenceCountWithoutBrazil->Draw();
-    c2->SaveAs((TString)outdir+"PIC003A_WithoutBrazilExpositionAndCount.jpg");
-    
-    c3->cd(1);
-    MainAxisPamelaExpositionEquator->Draw();
-    c3->cd(2);
-    MainAxisPamelaCountEquator->Draw();
-    c3->cd(3);
-    PitchExpositionNoOrientationEquator->Draw();
-    c3->cd(4);
-    PitchCountNoOrientationEquator->Draw();
-    c3->cd(5);
-    DiferenceBetweenNoOrientationAndWithOrientationEquator->Divide(PitchExpositionNoOrientationEquator,MainAxisPamelaExpositionEquator);
-    DiferenceBetweenNoOrientationAndWithOrientationEquator->Draw();
-    c3->cd(6);
-    for(Int_t iu=0;iu<360;iu++){
-	BinContent = DiferenceBetweenNoOrientationAndWithOrientationEquator->GetBinContent(iu);
-	DifferenceCountEquator->Fill(BinContent);
-    }
-    DifferenceCountEquator->Draw();
-    c3->SaveAs((TString)outdir+"PIC005A_EquatorExpositionAndCount.jpg");
-    
-    c4->cd(1);
-    PitchExpositionBrazil->Draw();
-    c4->cd(2);
-    PitchCountBrazil->Draw();
-    c4->cd(3);
-    PitchExpositionWithoutBrazil->Draw();
-    c4->cd(4);
-    PitchCountWithoutBrazil->Draw();
-    c4->cd(5);
-    PitchExpositionEquator->Draw();
-    c4->cd(6);
-    PitchCountEquator->Draw();
-    c4->cd(7);
-    ZenitAngle->Draw();
-    c4->SaveAs((TString)outdir+"PIC006A_PitchCountsAndExpositions.jpg");
-
-    c5->cd(1);
-    PitchvsLatInBrazil->Draw();
-    c5->cd(2);
-    LieveTimevsLatInBrazil->Draw();
-    c5->SaveAs((TString)outdir+"PIC006A_PitchvsLotitudeBrazil.jpg");
-
+    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();
+    	    }
+    	    
+        }
     }
-
-    /***************************************************************************************/
-    //Draw hystogramms for diferences in results between calculating using DoTrack2
-    //and using axv[0], ayv[0]
-    /**************************************************************************************/
-
-    if(DiffHist){
-
-    c6->cd(1);
-    DifNicoANdMeAzim->Draw();
-    c6->cd(2);
-    DifNicoANdMeZenit->Draw();
-    c6->SaveAs((TString)outdir+"PIC001Diference.jpg");
-    
+    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 hystogram to count how much time Pamela have each Pitch-Angle
-    /***********************************************************************/
-    
-    if(PamAngTime) {
-
-    c8->cd();
-    PitchTime->Draw();
-    c8->SaveAs((TString)outdir+"PIC002PamAngTime.jpg");
 
-    }
-    
-    if(PamExp){
-    
-    c9->cd();
-    PitchExposure->Draw();
-    c9->SaveAs((TString)outdir+"PIC003PamAngTime.jpg");
-    
-    }
-    
-    /**************************************************************************************/
-    //Draw Hystogramm for Pamela's angular efficiency
-    /**************************************************************************************/
-    gStyle->SetPalette(1);
-    
-    if(PamEff){
-    c7->cd(1);
-    PamAngEffhigh->Draw("colz");
-    TFile fh((TString)outdir+"PamAngEfficiencyHighEnergy.root","recreate");
-    PamAngEffhigh->Write();
-    fh.Close();
-    //PamAngEffhigh->SaveAs((TString)outdir+"PamAngEfficiencyHighEnergy.root");
-    c7->cd(2);
-    PamAngEfflow->Draw("colz");
-    TFile fl((TString)outdir+"PamAngEfficiencyLowEnergy.root","recreate");
-    PamAngEffhigh->Write();
-    fl.Close();
-    //PamAngEfflow->SaveAs((TString)outdir+"PamAngEfficiencyLowEnergy.root");
-    c7->SaveAs((TString)outdir+"PIC001PamAngEfficiencyHighEnergy.jpg");
-    
-    }
-    
     /**************************************************************************************/
     //Draw hystogramms for protons. Volodia
     /**************************************************************************************/