/[PAMELA software]/ParDirCalc/QtoInclination.C

Diff of /ParDirCalc/QtoInclination.C

Parent Directory | Revision Log | View Patch Patch

-revision 1.1.1.1 by cafagna,
Thu Jan 17 16:15:13 2008 UTC
+revision 1.3 by pam-mep,
Tue Sep 16 09:38:49 2008 UTC
 Line 6
  //
  #include "TrAng.h"
  #include "QtoInc.h"
+ #include "option.h"
  #include <fstream>
  #include <iostream>
-Line 24
+Line 25
  #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;
-Line 41 
 Bool_t DoTr = false;
+Line 397 
 Bool_t DoTr = false;
  Bool_t PamAngTime = false;
  Bool_t PamExp = false;
+ /*vector<TH1F> Hyst(1);
+ Hyst[0].SetName("ElectronLowEnergy");
+ Hyst[0].SetBins(36,0,360);
+ Hyst[0].FillRandom("gaus",5000);
+ TFile f1("/home/pamelaprod/malakhov/QtoInc/TEST/TEST.root","recreate");
+ Hyst[0].Write();
+ f1.Close();
+ //TFile* f1 = new TFile("/home/pamelaprod/malakhov/QtoInc/TEST/TEST.root");
+ //TH1F* h1 = (TH1F*)f1->Get("");
+ //cout<<h1->GetEntries()<<endl;
+ */
  if(argc<2){
      cout<<"You have to insert at least a file to analisys \nUsing --help for more information\n";
      exit(1);
  }
+ //opt1.opt.PE.Proton.Hyst1.resize(opt1.opt.PE.Proton.Hyst1.size()+1);
+ //opt1.opt.PE.Proton.Hyst1[0].SetBins(36,0,360);
+ //opt1.opt.PE.Proton.Hyst1[0].SetName("Test2");
+ //cout<<opt1.opt.PE.Proton.Hyst1.size()<<endl;
+ //cout<<opt1.opt.PE.Proton.Hyst1[0].GetName()<<endl;
+ //vector<TH2F>Hyst2DF(1);
+ //Hyst2DF[0].SetName("Test1");
+ //Hyst2DF[0].SetBins(36,0,360,36,0,360);
+ //Hyst2DF.resize(Hyst2DF.size()+1);
+ //Hyst2DF[1].SetName("Test2");
+ //Hyst2DF[1].SetBins(36,0,360,36,0,360);
+ //cout<<Hyst2DF[0].GetName()<<endl;
  if(!strcmp(argv[1],"--help")){
-     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";
-Line 61 
 if(!strcmp(argv[1],"--help")){
+Line 443 
 if(!strcmp(argv[1],"--help")){
      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(!strcmp(argv[1],"--Divide")){
-     if (4 >= argc){
+     if(argc<7){
-     cout<<"--OutPath needs argument\n See --help\n";
+         cout<<"--Divide  needs argument\n See --help\n";
-     }else outdir = argv[4];
+         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[1],"--ReBin")){
-     if(!strcmp(argv[i],"-DoTr")) DoTr = true;
+     if(4 >= argc){
-     if(!strcmp(argv[i],"-PamEff")) PamEff = true;
+         cout<<"--ReBin  needs argument\n See --help\n";
-     if(!strcmp(argv[i],"-DiffHist")){DiffHist = true; DoTr = true;}
+         exit(1);
-     if(!strcmp(argv[i],"-AngHist")) AnglHist = true;
+     }
-     if(!strcmp(argv[i],"-PamAngTime")) PamAngTime = true;
+     ReBinTH2F((TString)argv[2],(TString)argv[3],argv[4]);
-     if(!strcmp(argv[i],"-PamExp")) PamExp = true;
+     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
-Line 116 
 Float_t ConvertT2R(Float_t T, Float_t M,
+Line 1055 
 Float_t ConvertT2R(Float_t T, Float_t M,
   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);
+     if(opt.opt.verbose>=2) cout<<"\nGeneral variables initialisation...\n\n";
-     TSystemDirectory *workdir = new TSystemDirectory("workdir",dirname);
-     TList *flist=workdir->GetListOfFiles();
+     TSystemDirectory* workdir = new TSystemDirectory("workdir",dirname);
+     //TSystemDirectory* workdir2 = new TSystemDirectory("workdir2","/data01/_3/607_3/");
+     TList* flist = workdir->GetListOfFiles();
      PamLevel2* pam_events = new PamLevel2();
      PamelaOrientation* PO = new PamelaOrientation();
      TTree *T = pam_events->GetPamTree(flist,"treename");
      ULong_t nevents = T->GetEntries();
+     if(opt.opt.verbose>=1) cout<<"\n Nevents = "<<nevents<<"\n\n";
-     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;
+     Int_t ntr;                  Float_t Argv = 0;       Double_t PamAzim = 0;
-     Float_t Argv = 0;
+     Double_t PamZenit = 0;      Double_t MyAzim = 0;    Double_t MyZenit = 0;
-     Double_t PamAzim = 0;
+     Double_t Px = 0;            Double_t Py = 0;        Double_t Pz = 0;
-     Double_t PamZenit = 0;
-     Double_t MyAzim = 0;
+     //Peremennye dlya opredeleniya pervogo vhoda v cycly v PamExp
-     Double_t MyZenit = 0;
-     Double_t Px;
+     ULong_t PHE1 = 0;           ULong_t PLE1 = 0;       ULong_t PAE1 = 0;
-     Double_t Py;
+     ULong_t EHE1 = 0;           ULong_t ELE1 = 0;       ULong_t EAE1 = 0;
-     Double_t Pz;
+     ULong_t pHE1 = 0;           ULong_t pLE1 = 0;       ULong_t pAE1 = 0;
+     ULong_t AHE1 = 0;           ULong_t ALE1 = 0;       ULong_t AAE1 = 0;
-     /*********************************************************************************/
-     // 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");
-         }
-     }
      /*********************************************************************************/
      //Histogramms for Count and exposition of Angles (Pitch, Pamela's Main axis, etc )
      /*********************************************************************************/
-     //if(AnglHist){
+     if(opt.opt.verbose>=2) cout<<"\n--AnglHyst Hystogramm block...\n\n";
-     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);
-     //}
      /***************************************************************************************/
-     //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){
+     if(opt.opt.verbose>=2) cout<<"\n--DiffHyst Hystogramm block...\n\n";
-     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);
      /***************************************************************************************/
      //Initialization histogramms for count of protons depends from Energy. Volodia.
      /***************************************************************************************/
- /*
-     TH1F* proton10log= new TH1F("Protons Flux","",80,-1.,3.);
+ //    TH1F* proton10log= new TH1F("Protons Flux","",80,-1.,3.);
-     proton10log->GetXaxis()->SetTitle("log10(Energy) , GeV");
+ //    proton10log->GetXaxis()->SetTitle("log10(Energy) , GeV");
-     TH1F* proton1log= new TH1F("Protons Flux","",80,-1.,3.);
+ //    TH1F* proton1log= new TH1F("Protons Flux","",80,-1.,3.);
-     TH1F* binw= new TH1F("w","",80,-1.,3.);
+     //TH1F* binw= new TH1F("w","",80,-1.,3.);
-     proton1log->GetXaxis()->SetTitle("log10(Energy) , GeV");
+     TH1F binw;
-     proton1log->SetLineColor(kRed);
+     binw.SetBins(80,-1.,3.);
-     TCanvas *plogCanvas = new TCanvas("proton flux","protonflux", 800, 800);
+ //    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 x;
-     Float_t xmin, xmax;
+     //Float_t xmin, xmax;
-     Float_t binwide[80];
+     //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++) {
+     //for (Int_t l=0 ; l<80; l++) {
-     binwide[l]=binw->GetBinContent(l+1);
+     //binwide[l]=binw->GetBinContent(l+1);
-     }
+ //    }
- */
      /****************************************************************************/
      //Geometrical Factor. Volodia
      /****************************************************************************/
-Line 305 
 void Calculate(char* dirname, string out
+Line 1161 
 void Calculate(char* dirname, string out
      //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){
+         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
-             //Getting general parameters
+         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
-             track = pam_events->GetTrack(0);
+         Float_t Alt = pam_events->GetOrbitalInfo()->alt;
-             Int_t M0DE = pam_events->GetOrbitalInfo()->mode;  //0 is zero here
+         Float_t Lat = pam_events->GetOrbitalInfo()->lat;
-             Float_t Bx = -pam_events->GetOrbitalInfo()->Bdown;
+         Float_t Lon = pam_events->GetOrbitalInfo()->lon;
-             Float_t By = pam_events->GetOrbitalInfo()->Beast;
+         Float_t LieveTime = 0.16*pam_events->GetTrigLevel2()->dltime[0];
-             Float_t Bz = pam_events->GetOrbitalInfo()->Bnorth;
+         Float_t DeathTime = 0.01*pam_events->GetTrigLevel2()->dltime[1];
-             Float_t Babs = pam_events->GetOrbitalInfo()->Babs;
+         Float_t rigev = 0;
-             Float_t L = pam_events->GetOrbitalInfo()->L;
+         Float_t detr = 0;
-             Float_t rigev = 1/track->GetTrkTrack()->GetDeflection();
+         Int_t hm0; Int_t hm1;
+         Int_t hs0; Int_t hs1;
-             /*************************************************************************/
-             //Track selection
+         ntr = 0;
-             /*************************************************************************/
+         Double_t A1; Double_t A2; Double_t A3;
-             ntr = 0;
+         Double_t SB = 1000; Double_t SA; Double_t SC; Double_t ZC;
-             if ((track->GetTrkTrack()->chi2<=0) || (track->GetTrkTrack()->chi2>=500)) ntr=-1;
-             for(Int_t ii=1;ii<=12;ii++){
+         if(OneTrack){
-                 if ((track->GetToFTrack()->beta[i]<0) || (track->GetToFTrack()->beta[i]>99)) ntr=-1;
+             track = pam_events->GetTrack(0);                            //don't need for PamExp option
-             }
+             rigev = 1/track->GetTrkTrack()->GetDeflection();            //dont need for AnglHyst, DiffHyst and PamExp options;
-             if((track->GetTrkTrack()->GetNX()<=4)&&(track->GetTrkTrack()->GetNY()<4)) ntr=-1;
+             detr = track->GetTrkTrack()->GetDEDX();
-             if((M0DE!=0)&&(M0DE!=1)&&(M0DE!=6)&&(M0DE!=2)&&(M0DE!=3)&&(M0DE!=8)&&(M0DE!=4)) ntr=-1;
+             hm0 = pam_events->GetAcLevel2()->hitmap[0];
-             if(sqrt(pow(Bx,2)+pow(By,2)+pow(Bz,2))<0.001) ntr=-1;
+             hm1 = pam_events->GetAcLevel2()->hitmap[1];
-             if(rigev<0.600) ntr=-1;
+             hs0 = pam_events->GetAcLevel2()->hitstatus[0];
+             hs1 = pam_events->GetAcLevel2()->hitstatus[1];
-             if(ntr!=-1 || PamExp){
+             if(!OneTrack) ntr = -1; else{
-                 /*******************************************************************************/
+                 //cout<<"OneTrack"<<endl;
-                 //Transit reference system
+                 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++){
-                 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);
+                     if ((track->GetToFTrack()->beta[ii]<0) || (track->GetToFTrack()->beta[ii]>99)) ntr=-1;
-                 TMatrixD Dij = PO->GreenwichtoGEO(pam_events->GetOrbitalInfo()->lat,pam_events->GetOrbitalInfo()->lon,Fij);
+                 }*/
-                 TMatrixD Iij = PO->ColPermutation(Dij);
+                 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;
-                 //Calculate Zenit and Azimutal angle and vector of particle flight in Pamela
+                 if(sqrt(pow(Bx,2)+pow(By,2)+pow(Bz,2))<0.001) ntr=-1;
-                 //using axv[0] and ayv[0] variables
+                 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();
                  PamZenit = TMath::RadToDeg()*asin(sqrt(pow(sin(Aayv), 2) + pow(sin(Aaxv), 2)));
                  Double_t axv = -track->GetTrkTrack()->axv[0] * TMath::DegToRad();
                  Double_t ayv = -track->GetTrkTrack()->ayv[0] * TMath::DegToRad();
                  Double_t angle = atan(sin(TMath::Abs(ayv))/sin(TMath::Abs(axv))) * TMath::RadToDeg();
                  PamAzim =  360. - angle;
                  if(axv>=0 && ayv >=0) PamAzim = angle;
                  if(axv<0 && ayv >0) PamAzim = 180. - angle;
                  if(axv<0 && ayv <0) PamAzim = 180. + angle;
                  PamAzim = PamAzim * TMath::DegToRad();
                  PamZenit = (180 - PamZenit) * TMath::DegToRad();
                  Px = sin(axv);
                  Py = sin(ayv);
                  Pz = cos(PamZenit);
+                 //cout<<"PamAzimuth = "<<PamAzim<<" PamZenit = "<<PamZenit<<endl;
-                 }
+             }
-                 /*****************************************************************************/
+             /*****************************************************************************/
-                 //Calculate Zenit and Azimutal angle and vector of particle flight in Pamela
+             //Calculate Zenit and Azimutal angle and vector of particle flight in Pamela
-                 //using DoTrack2 procedure
+             //using DoTrack2 procedure
-                 /*****************************************************************************/
+             /*****************************************************************************/
-                 if(DoTr){
+             if(opt.opt.DoTr){
                  track->GetTrkTrack()->DoTrack2(tr);
                  Double_t E11x = tr->x[0];
                  Double_t E11y = tr->y[0];
                  Double_t E11z = zin[0];
                  Double_t E22x = tr->x[3];
                  Double_t E22y = tr->y[3];
                  Double_t E22z = zin[3];
                  Double_t norm = sqrt(pow(E22x-E11x,2)+pow(E22y-E11y,2)+pow(E22z-E11z,2));
                  MyAzim = TMath::RadToDeg()*atan(TMath::Abs(E22y-E11y)/TMath::Abs(E22x-E11x));
                  MyAzim = MyAzim;//-180;
                  MyZenit = TMath::RadToDeg()*acos((E22z-E11z)/norm);
                  if(E22x-E11x>=0 && E22y-E11y <0) MyAzim =  360. - MyAzim;
                  if(E22x-E11x>=0 && E22y-E11y >=0) MyAzim = MyAzim;
                  if(E22x-E11x<0 && E22y-E11y >0) MyAzim = 180. - MyAzim;
                  if(E22x-E11x<0 && E22y-E11y <0) MyAzim = 180. + MyAzim;
                  Px = (E22x-E11x)/norm;
                  Py = (E22y-E11y)/norm;
                  Pz = (E22z-E11z)/norm;
-                 }
-                 /***************************************************************************************/
-                 //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.
-                 /******************************************************************************/
+             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);
-                 if(14.9/L/L > 0. && 14.9/L/L < 1. & Babs >0.24) {
+             TMatrixD Iij = PO->ColPermutation(Dij);
-                 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){
-                 if(L>6){
+             TMatrixD Eij = PO->PamelatoGEO(Iij,Px,Py,Pz);               //Don't necessary for PamExp
-                 Float_t Mass = 0;
-                 Int_t detr = 0;
+             A1 = Iij(0,2);
-                 Float_t Ekin = 0;
+             A2 = Iij(1,2);
-                 Float_t betaev = 0.25*(track->GetToFTrack()->beta[0]+track->GetToFTrack()->beta[1]+track->GetToFTrack()->beta[2]+track->GetToFTrack()->beta[3]);
+             A3 = Iij(2,2);
-                     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) ) {
+             SB = PO->GetPitchAngle(Eij(0,0),Eij(1,0),Eij(2,0),Bx,By,Bz);        // Pitch angle
-                             Ekin=ConvertR2T(rigev,0.511,-1.);
+             SA = PO->GetPitchAngle(A1,A2,A3,Bx,By,Bz);                  // Angle between Pamela's main axiz and B
-                             if (Ekin>1) PamAngEffhigh->Fill(TMath::RadToDeg()*PamAzim,TMath::RadToDeg()*PamZenit);
+             SC = PO->GetPitchAngle(1,0,0,Bx,By,Bz);                             // Angle between direction to zenit and B
-                             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.);
+         //cout<<"OneTrack"<<endl;
-                             if (Ekin>1) PamAngEffhigh->Fill(TMath::RadToDeg()*PamAzim,TMath::RadToDeg()*PamZenit);
+         if(sqrt(pow(Bx,2)+pow(By,2)+pow(Bz,2))<0.001) ntr=-1;else{
-                             if (Ekin<1) PamAngEfflow->Fill(TMath::RadToDeg()*PamAzim,TMath::RadToDeg()*PamZenit);
+             //cout<<"US.size = "<<opt.US.size()<<endl;
+             for(Int_t jz = 0; jz<opt.PS.size();jz++){
+                 for(Int_t jy = 0; jy<opt.US.size();jy++){
+                     for(Int_t jx = 0; jx<opt.US[jy].USP.size();jx++){
+                         //if(SB<400) cout<<"Pitch = "<<SB<<" Babs = "<<Babs<<" L = "<<L<<" Alt = "<<Alt<<endl;
+                         //cout<<"Fill Hyst1D["<<thg<<"]"<<endl;
+                         if(Babs<opt.US[jy].USP[jx].Bless && Babs>opt.US[jy].USP[jx].Bmore && L<opt.US[jy].USP[jx].Lless && L>opt.US[jy].USP[jx].Lmore && Alt<opt.US[jy].USP[jx].altless && Alt>opt.US[jy].USP[jx].altmore  && SB<opt.PS[jz].PitchMax && SB>opt.PS[jz].PitchMin){
+                             for(Int_t hg = 0; hg<opt.US[jy].NAH1D.size(); hg++){
+                                 Int_t thg = opt.US[jy].NAH1D[hg];
+                                 for(Int_t ps = 0; ps<opt.PS[jz].NAH1D.size();ps++){
+                                     if(thg==opt.PS[jz].NAH1D[ps]){
+                                         opt.Hyst1D[thg].HI.LieveTime+= 0.16*pam_events->GetTrigLevel2()->dltime[0];
+                                         opt.Hyst1D[thg].HI.DeathTime+= 0.01*pam_events->GetTrigLevel2()->dltime[1];
+                                         //cout<<"Pitch = "<<SB<<" Babs = "<<Babs<<" L = "<<L<<" Alt = "<<Alt<<endl;
+                                         //cout<<"Fill Hyst1D["<<thg<<"]"<<endl;
+                                     }
+                                 }
+                             }
+                             for(Int_t hg = 0; hg<opt.US[jy].NAH2D.size(); hg++){
+                                 Int_t thg = opt.US[jy].NAH2D[hg];
+                                 for(Int_t ps = 0; ps<opt.PS[jz].NAH2D.size();ps++){
+                                     if(thg==opt.PS[jz].NAH2D[ps]){
+                                         opt.Hyst2D[thg].HI.LieveTime+= 0.16*pam_events->GetTrigLevel2()->dltime[0];
+                                         opt.Hyst2D[thg].HI.DeathTime+= 0.01*pam_events->GetTrigLevel2()->dltime[1];
+                                     }
+                                 }
+                             }
+                             for(Int_t hg = 0; hg<opt.US[jy].NAH3D.size(); hg++){
+                                 Int_t thg = opt.US[jy].NAH3D[hg];
+                                 for(Int_t ps = 0; ps<opt.PS[jz].NAH3D.size();ps++){
+                                     if(thg==opt.PS[jz].NAH3D[ps]){
+                                         opt.Hyst3D[thg].HI.LieveTime+= 0.16*pam_events->GetTrigLevel2()->dltime[0];
+                                         opt.Hyst3D[thg].HI.DeathTime+= 0.01*pam_events->GetTrigLevel2()->dltime[1];
+                                     }
+                                 }
+                             }
+                             jx = opt.US[jy].USP.size();
                          }
                      }
                  }
+             }
-                 }
+         }
+         //cout<<ntr<<endl;
-                 /***********************************************************************/
+         if(OneTrack){
-                 //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){
-                 /***********************************************************************/
+             if(ntr!=-1){
-                 //Filling angular histogramms for Brazil
-                 /***********************************************************************/
+                 Bool_t verb = false;
+                 if(opt.opt.Nverb!=-1) if(i>=opt.opt.Nverb) verb = true;
-                 if((Babs<0.19)&&(L<1.75)){
+                 //if(opt.opt.Nverb==-1) verb = false;
-                         Argv = 0.16*pam_events->GetTrigLevel2()->dltime[0];
-                         PitchExpositionBrazil->Fill(SB,Argv);
+                 if(verb) cout<<"\n\ni = "<<i<<"\n\n"<<endl;
-                         MainAxisPamelaExpositionBrazil->Fill(SA,Argv);
-                         PitchExpositionNoOrientationBrazil->Fill(SC,Argv);
+                 /*****************************************************************************/
-                         PitchCountBrazil->Fill(SB);
+                 // CALO Selection
-                         MainAxisPamelaCountBrazil->Fill(SA);
+                 /*****************************************************************************/
-                         PitchCountNoOrientationBrazil->Fill(SC);
-                         ZenitAngle->Fill(ZC);
+                 Float_t caloqpre=track->GetCaloTrack()->qpre;
-                         PitchvsLatInBrazil->Fill(pam_events->GetOrbitalInfo()->lat,SA);
+                 Int_t   calonpre=track->GetCaloTrack()->npre;
-                         LieveTimevsLatInBrazil->Fill(pam_events->GetOrbitalInfo()->lat,Argv);
+                 Float_t caloqtrack=track->GetCaloTrack()->qtrack;
-                 }
+                 Float_t caloqmax=pam_events->GetCaloLevel2()->qmax;
+                 Int_t   dncore=track->GetCaloTrack()->ncore;
-                 /***********************************************************************/
+                 Int_t   lepton = 0;
-                 //Filling angulars histogramms for areas without Brazil
+                 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){
+                 if(verb){
-                         Argv = 0.16*pam_events->GetTrigLevel2()->dltime[0];
+                 //if(rigev>1.7 && rigev<= 1.8 && L>6 && Lat>60 && Lat<-60){
-                         PitchExpositionWithoutBrazil->Fill(SB,Argv);
+                     cout<<"rigev = "<<rigev<<"\tdetr = "<<detr<<endl;
-                         MainAxisPamelaExpositionWithoutBrazil->Fill(SA,Argv);
+                     cout<<"Babs = "<<Babs<<"\tL = "<<L<<endl;
-                         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);
-                         }
                  }
+                 /*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)
+                 Double_t Rmin = 0;
+                 Double_t Rmax = 0;
-             } // if(PamEff||AnglHist)
+                 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;
-     /*********************************************************/
+     TCanvas* Canv;
-     //Divide canvases for angular histogramms
+     TCanvas* Canvdiv;
-     /*********************************************************/
+     TCanvas* CanvExp;
+     TCanvas* CanvLog;
-     if(AnglHist){
+     TString Ext[3];
+     Ext[0] = ".jpg";    Ext[1] = ".png";        Ext[2] = ".bmp";
-     c1->Divide(1,6);
-     c2->Divide(1,6);
-     c3->Divide(1,6);
-     c4->Divide(1,7);
-     c5->Divide(1,2);
-     }
+     if(opt.opt.verbose>=2) cout<<"\n\nSave block...\n\n";
-     /***************************************************************************************/
-     //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(PamEff){
+     gStyle->SetPalette(1);
-     c7->Divide(1,2);
-     }
-     /**************************************************************************************/
+     for(Int_t i = 0; i<opt.Hyst1D.size(); i++){
-     //Draw Angular Histogramms
+         TH1F hdiv;
-     /**************************************************************************************/
+         TH1F h1n1D; TH1F h2n1D; TH1F hlog;
+         if(opt.Hyst1D[i].HI.ExpType!=0){
-     if(AnglHist){
+             hlog.SetBins(80,-1.,3.);
+             h1n1D.SetBins(opt.Hyst1D[i].Hyst1DF.GetNbinsX(),opt.Hyst1D[i].Hyst1DF.GetXaxis()->GetXmin(),opt.Hyst1D[i].Hyst1DF.GetXaxis()->GetXmax());
-     c1->cd(1);
+             h2n1D.SetBins(opt.Hyst1D[i].Hyst1DF.GetNbinsX(),opt.Hyst1D[i].Hyst1DF.GetXaxis()->GetXmin(),opt.Hyst1D[i].Hyst1DF.GetXaxis()->GetXmax());
-     MainAxisPamelaExpositionBrazil->Draw();
+             //h2div.SetBins(h12D->GetNbinsX(),h12D->GetXaxis()->GetXmin(),h12D->GetXaxis()->GetXmax(),h12D->GetNbinsY(),h12D->GetYaxis()->GetXmin(),h12D->GetYaxis()->GetXmax());
-     c1->cd(2);
+             for(Int_t ii = 0; ii < opt.Hyst1D[i].Hyst1DF.GetNbinsX(); ii++){
-     MainAxisPamelaCountBrazil->Draw();
+                 h1n1D.SetBinContent(ii,opt.Hyst1D[i].Hyst1DF.GetBinContent(ii)/opt.Hyst1D[i].Hyst1DF.GetEntries());
-     c1->cd(3);
+                 h2n1D.SetBinContent(ii,opt.Hyst1D[i].Hyst1DFExp.GetBinContent(ii)/opt.Hyst1D[i].Hyst1DFExp.GetEntries());
-     PitchExpositionNoOrientationBrazil->Draw();
+             }
-     c1->cd(4);
+             hdiv.SetName(opt.Hyst1D[i].Hyst1DF.GetName());
-     PitchCountNoOrientationBrazil->Draw();
+             hdiv.SetBins(opt.Hyst1D[i].Hyst1DF.GetNbinsX(),opt.Hyst1D[i].Hyst1DF.GetXaxis()->GetXmin(),opt.Hyst1D[i].Hyst1DF.GetXaxis()->GetXmax());
-     c1->cd(5);
+             hdiv.Divide(&h1n1D,&h2n1D);
-     DiferenceBetweenNoOrientationAndWithOrientationBrazil->Divide(PitchExpositionNoOrientationBrazil,MainAxisPamelaExpositionBrazil);
+             cout<<"i = "<<i<<endl;
-     DiferenceBetweenNoOrientationAndWithOrientationBrazil->Draw();
+             if(opt.Hyst1D[i].XType==17) hlog.Divide(&opt.Hyst1D[i].Hyst1DF,&binw);
-     c1->cd(6);
+         }
-     for(Int_t iu=0;iu<360;iu++){
+         if(opt.Hyst1D[i].HI.TypeExt > 0){
-         BinContent = DiferenceBetweenNoOrientationAndWithOrientationBrazil->GetBinContent(iu);
+             Canv = new TCanvas("Canv",opt.Hyst1D[i].Hyst1DF.GetName(),1200,800);
-         DifferenceCountInBrazil->Fill(BinContent);
+             Canv->cd();
-     }
+             opt.Hyst1D[i].Hyst1DF.Draw();
-     DifferenceCountInBrazil->Draw();
+             Canv->SaveAs((TString)outdir+"ParDirCalc_"+opt.Hyst1D[i].HI.filename+".png");
-     //DifferenceCountInBrazil->SaveAs((TString)outdir+"PIC001A_BrazilExpositionAndCount.root");
+             if(opt.Hyst1D[i].HI.ExpType!=0){
-     c1->SaveAs((TString)outdir+"PIC001A_BrazilExpositionAndCount.jpg");
+                 CanvExp = new TCanvas("CanvExp",opt.Hyst1D[i].Hyst1DFExp.GetName(),1200,800);
+                 CanvExp->cd();
-     c2->cd(1);
+                 opt.Hyst1D[i].Hyst1DFExp.Draw();
-     MainAxisPamelaExpositionWithoutBrazil->Draw();
+                 CanvExp->SaveAs((TString)outdir+"ParDirCalc_"+opt.Hyst1D[i].HI.filename+"_Exp.png");
-     c2->cd(2);
+                 Canvdiv = new TCanvas("Canvdiv",hdiv.GetName(),1200,800);
-     MainAxisPamelaCountWithoutBrazil->Draw();
+                 Canvdiv->cd();
-     c2->cd(3);
+                 hdiv.Draw();
-     PitchExpositionNoOrientationWithoutBrazil->Draw();
+                 Canvdiv->SaveAs((TString)outdir+"ParDirCalc_"+opt.Hyst1D[i].HI.filename+"DivExp.png");
-     c2->cd(4);
+             }
-     PitchCountNoOrientationWithoutBrazil->Draw();
+             if(opt.Hyst1D[i].XType==17){
-     c2->cd(5);
+                 CanvLog = new TCanvas("CanvLog",opt.Hyst1D[i].Hyst1DF.GetName(),1200,800);
-     DiferenceBetweenNoOrientationAndWithOrientationWithoutBrazil->Divide(PitchExpositionNoOrientationWithoutBrazil,MainAxisPamelaExpositionWithoutBrazil);
+                 CanvLog->cd();
-     DiferenceBetweenNoOrientationAndWithOrientationWithoutBrazil->Draw();
+                 hlog.Draw();
-     c2->cd(6);
+                 CanvLog->SaveAs((TString)outdir+"ParDirCalc_"+opt.Hyst1D[i].HI.filename+"_DivBinWide.png");
-     for(Int_t iu=0;iu<360;iu++){
+             }
-         BinContent = DiferenceBetweenNoOrientationAndWithOrientationWithoutBrazil->GetBinContent(iu);
+         }
-         DifferenceCountWithoutBrazil->Fill(BinContent);
+         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");
-     DifferenceCountWithoutBrazil->Draw();
+             opt.Hyst1D[i].Hyst1DF.Write();
-     c2->SaveAs((TString)outdir+"PIC003A_WithoutBrazilExpositionAndCount.jpg");
+             fg.Close();
+             TFile fe((TString)outdir+"ParDirCalc_"+opt.Hyst1D[i].HI.filename+"_DivWide.root","recreate");
-     c3->cd(1);
+             hlog.Write();
-     MainAxisPamelaExpositionEquator->Draw();
+             fe.Close();
-     c3->cd(2);
+             if(opt.Hyst1D[i].HI.ExpType!=0){
-     MainAxisPamelaCountEquator->Draw();
+                 TFile fgExp((TString)outdir+"ParDirCalc_"+opt.Hyst1D[i].HI.filename+"_Exp.root","recreate");
-     c3->cd(3);
+                 opt.Hyst1D[i].Hyst1DFExp.Write();
-     PitchExpositionNoOrientationEquator->Draw();
+                 fg.Close();
-     c3->cd(4);
+                 TFile fgdiv((TString)outdir+"ParDirCalc_"+opt.Hyst1D[i].HI.filename+"_DivExp.root","recreate");
-     PitchCountNoOrientationEquator->Draw();
+                 hdiv.Write();
-     c3->cd(5);
+                 fgdiv.Close();
-     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.Hyst2D.size(); i++){
-     /***************************************************************************************/
+         TH2F hdiv;
-     //Draw hystogramms for diferences in results between calculating using DoTrack2
+         TH2F h1n2D; TH2F h2n2D;
-     //and using axv[0], ayv[0]
+         cout<<"i = "<<i<<endl;
-     /**************************************************************************************/
+         if(opt.Hyst2D[i].HI.ExpType!=0){
+             //cout<<"Set hdiv"<<endl;
-     if(DiffHist){
+             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());
-     c6->cd(1);
+             //h2div.SetBins(h12D->GetNbinsX(),h12D->GetXaxis()->GetXmin(),h12D->GetXaxis()->GetXmax(),h12D->GetNbinsY(),h12D->GetYaxis()->GetXmin(),h12D->GetYaxis()->GetXmax());
-     DifNicoANdMeAzim->Draw();
+             for(Int_t ii = 0; ii < opt.Hyst2D[i].Hyst2DF.GetNbinsX(); ii++){
-     c6->cd(2);
+                 for(Int_t j = 0; j < opt.Hyst2D[i].Hyst2DF.GetNbinsY(); j++){
-     DifNicoANdMeZenit->Draw();
+                     h1n2D.SetBinContent(ii,j,opt.Hyst2D[i].Hyst2DF.GetBinContent(ii,j)/opt.Hyst2D[i].Hyst2DF.GetEntries());
-     c6->SaveAs((TString)outdir+"PIC001Diference.jpg");
+                     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
      /**************************************************************************************/

 Legend:



Removed from v.1.1.1.1
 


changed lines


 
Added in v.1.3
 Legend:



Removed from v.1.1.1.1
 


changed lines


 
Added in v.1.3
-Removed from v.1.1.1.1
+Added in v.1.3

	ViewVC Help
Powered by ViewVC 1.1.23