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

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

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

using namespace std;
     
int main(int argc, char* argv[]){

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

if(argc<2){
    cout<<"You have to insert at least a file to analisys \nUsing --help for more information\n";
    exit(1);
}

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

if(!strcmp(argv[3],"--OutPuth")){
    if (4 >= argc){
    cout<<"--OutPath needs argument\n See --help\n";    
    }else outdir = argv[4];
}

for(Int_t i=4; i<argc; i++){
    if(!strcmp(argv[i],"-DoTr")) DoTr = true;
    if(!strcmp(argv[i],"-PamEff")) PamEff = true;
    if(!strcmp(argv[i],"-DiffHist")){DiffHist = true; DoTr = true;}
    if(!strcmp(argv[i],"-AngHist")) AnglHist = true;
    if(!strcmp(argv[i],"-PamAngTime")) PamAngTime = true;
    if(!strcmp(argv[i],"-PamExp")) PamExp = true;
}

if(!PamEff && !DiffHist && !AnglHist && !PamAngTime && !PamExp) AnglHist = true;

if(!strcmp(argv[1],"--WorkDir")){
    if (2 >= argc){
    cout<<"--WorkDir needs argument\n See --help\n";
    exit(1);
    }else Calculate(argv[2],outdir,AnglHist,DiffHist,PamEff,DoTr,PamAngTime,PamExp);
}

}

//Function to Convert Rigidity to Energy

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

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

//Function to convert Energy to Rigidity

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

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

void Calculate(char* dirname, string outdir, Bool_t AnglHist, Bool_t DiffHist, Bool_t PamEff, Bool_t DoTr, Bool_t PamAngTime, Bool_t PamExp){
    
    /**********************************************/
    //First initialization, general for all purposes
    /**********************************************/
    
    //gStyle->SetOptStat(111111);
    TSystemDirectory *workdir = new TSystemDirectory("workdir",dirname);
    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();
    
    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/");
    /********************************************************************************/
    
    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);
    
    PamTrack *track;
    
    Int_t ntr;
    Float_t Argv = 0;
    Double_t PamAzim = 0;
    Double_t PamZenit = 0;
    Double_t MyAzim = 0;
    Double_t MyZenit = 0;
    Double_t Px;
    Double_t Py;
    Double_t Pz;
    
    /*********************************************************************************/
    // Angular exposure. How much time Pamela have any Azimutal and Zenital angle
    /*********************************************************************************/
    
    TH2F *hhigh;
    TH2F *hlow;
    
    TFile fhigh((TString)outdir+"PamAngEfficiencyHighEnergy.root");
    TFile flow((TString)outdir+"PamAngEfficiencyLowEnergy.root");
    if(PamExp){
    
	if (fhigh.IsZombie()||flow.IsZombie()){
	    cout<<"Problem with Hystogrammfiles"<<endl;
	    exit(1);
	}else{
	    hhigh = (TH2F*)fhigh.Get("PamAngEffhigh");
	    hlow = (TH2F*)flow.Get("PamAngEfflow");
	}
    
    }
    
    /*********************************************************************************/
    //Histogramms for Count and exposition of Angles (Pitch, Pamela's Main axis, etc )
    /*********************************************************************************/
    
    //if(AnglHist){
    
    TH1F *PitchExpositionWithoutBrazil = new TH1F("PitchExpositionWithoutBrazil", "Pitch Exposition without Brazil with Track", 360, 0, 180);
    TH1F *PitchExpositionEquator = new TH1F("PitchExpositionEquator", "Pitch Exposition in Equator  with Track", 360, 0, 180);
    TH1F *PitchExpositionBrazil = new TH1F("PitchExpositionBrazil", "Pitch Exposition Brazil with Track", 360, 0, 180);
    TH1F *MainAxisPamelaExpositionWithoutBrazil = new TH1F("MainAxisPamelaPitchExpositionWithoutBrazil", "Main Axis of Pamela Exposition without Brazil with Track", 360, 0, 180);
    TH1F *MainAxisPamelaExpositionEquator = new TH1F("MainAxisPamelaPitchExpositionEquator", "Main Axis of Pamela Exposition in Equator with Track", 360, 0, 180);
    TH1F *MainAxisPamelaExpositionBrazil = new TH1F("MainAxisPamelaExpositionBrazil", "Main Axis of Pamela Exposition Brazil with Track", 360, 0, 180);
    TH1F *PitchExpositionNoOrientationWithoutBrazil = new TH1F("PitchExpositionNoOrientationWithoutBrazil", "Pitch Exposition without orientation Information without Brazil with Track", 360, 0, 180);
    TH1F *PitchExpositionNoOrientationEquator = new TH1F("PitchExpositionNoOrientationEquator", "Pitch Exposition without orientation Information in Equator with Track", 360, 0, 180);
    TH1F *PitchExpositionNoOrientationBrazil = new TH1F("PitchExpositionNoOrientationBrazil", "Pitch Exposition without orientation Information in Brazil with Track", 360, 0, 180);
    TH1F *DiferenceBetweenNoOrientationAndWithOrientationWithoutBrazil = new TH1F("DiferenceBetweenNoOrientationAndWithOrientationWithoutBrazil", "Diference between no orientation data and with orientation without Brazil with Track", 360, 0, 180);
    TH1F *DiferenceBetweenNoOrientationAndWithOrientationEquator = new TH1F("DiferenceBetweenNoOrientationAndWithOrientationEquator", "Diference between no orientation data and with orientation in Equator with Track", 360, 0, 180);
    TH1F *DiferenceBetweenNoOrientationAndWithOrientationBrazil = new TH1F("DiferenceBetweenNoOrientationAndWithOrientationBrazil", "Diference between no orientation data and with orientation in Brazil with Track", 360, 0, 180);
    TH1F *DifferenceCountWithoutBrazil = new TH1F("DiferenceCountWithoutBrazil","Diference without Brazil",150,0.001,3);
    TH1F *DifferenceCountEquator = new TH1F("DiferenceCountWithoutBrazil","Diference without Brazil",150,0.001,3);
    TH1F *DifferenceCountInBrazil = new TH1F("DiferenceCountInBrazil","Diference in Brazil",150,0.001,3);
    
    TH1F *PitchCountWithoutBrazil = new TH1F("PitchCountWithoutBrazil", "Pitch Count without Brazil with Track", 360, 0, 180);
    TH1F *PitchCountEquator = new TH1F("PitchCountEquator", "Pitch Count in Equator with Track", 360, 0, 180);
    TH1F *PitchCountBrazil = new TH1F("PitchCountBrazil", "Pitch Count Brazil with Track", 360, 0, 180);
    TH1F *MainAxisPamelaCountWithoutBrazil = new TH1F("MainAxisPamelaPitchCountWithoutBrazil", "Main Axis of Pamela Count without Brazil with Track", 360, 0, 180);
    TH1F *MainAxisPamelaCountEquator = new TH1F("MainAxisPamelaPitchCountEquator", "Main Axis of Pamela Count in Equator with Track", 360, 0, 180);
    TH1F *MainAxisPamelaCountBrazil = new TH1F("MainAxisPamelaCountBrazil", "Main Axis of Pamela Count Brazil with Track", 360, 0, 180);
    TH1F *PitchCountNoOrientationWithoutBrazil = new TH1F("PitchCountNoOrientationWithoutBrazil", "Pitch Count without orientation Information without Brazil with Track", 360, 0, 180);
    TH1F *PitchCountNoOrientationEquator = new TH1F("PitchCountNoOrientationEquator", "Pitch Count without orientation Information in Equator with Track", 360, 0, 180);
    TH1F *PitchCountNoOrientationBrazil = new TH1F("PitchCountNoOrientationBrazil", "Pitch Count without orientation Information in Brazil with Track", 360, 0, 180);
    //TH1F *Diference2BetweenNoOrientationAndWithOrientationWithoutBrazil = new TH1F("Diference2BetweenNoOrientationAndWithOrientationWithoutBrazil", "Diference between no orientation data and with orientation without Brazil with Track", 360, 0, 180);
    //TH1F *Diference2BetweenNoOrientationAndWithOrientationEquator = new TH1F("Diference2BetweenNoOrientationAndWithOrientationEquator", "Diference between no orientation data and with orientation in Equator with Track", 360, 0, 180);
    //TH1F *Diference2BetweenNoOrientationAndWithOrientationBrazil = new TH1F("Diference2BetweenNoOrientationAndWithOrientationBrazil", "Diference between no orientation data and with orientation in Brazil with Track", 360, 0, 180);
    //TH1F *Difference2CountWithoutBrazil = new TH1F("DiferenceCountWithoutBrazil","Diference without Brazil",150,0.001,3);
    //TH1F *Difference2CountEquator = new TH1F("DiferenceCountWithoutBrazil","Diference without Brazil",150,0.001,3);
    //TH1F *Difference2CountInBrazil = new TH1F("DiferenceCountInBrazil","Diference in Brazil",150,0.001,3);
    //TH1F *EastPitchAllCount = new TH1F("EastPitchAllCount","East Count",)
    
    TH1F *ZenitAngle = new TH1F("ZenitAngle","Zenit Angle",360,120,180);
    //TH1F *ZenitAngleEquator = new TH1F("ZenitAngleWithoutEquator","Zenit Angle in Equator",360,120,180);
    //TH1F *ZenitAngleBrazil = new TH1F("ZenitAngleBrazil","Zenit Angle in Brazil",360,120,180);
    
    TH2F *PitchvsLatInBrazil = new TH2F("PitchvsLatInBrazil","Pitch Angle, deg; Latitude, deg",360,-40,-5,360,40,110);
    TH2F *LieveTimevsLatInBrazil = new TH2F("LieveTimevsLatInBrazil","Live time, deg; Latitude deg",360,-40,-5,1000,0,80);
    
    
    TCanvas *c1 = new TCanvas("c1","Exposition and count in Brazil",1200,3000);
    TCanvas *c2 = new TCanvas("c2","Exposition and Count without Brazil",1200,3000);
    TCanvas *c3 = new TCanvas("c3","Exposition and Count on Equator",1200,3000);
    TCanvas *c4 = new TCanvas("c4","Pitch-Angle exposition and count",1200,3400);
    TCanvas *c5 = new TCanvas("c5","PitchAndLieveTimevsLatitude",1200,800);
    TCanvas *c6 = new TCanvas("c6","Diferense",1200,800);
    //TCanvas *c6 = new TCanvas("c6","CountEquator",1200,3000);
    
    //}
    
    /***************************************************************************************/
    //Histogramms for diferences in results between calculating using DoTrack2
    //and using axv[0], ayv[0]
    /**************************************************************************************/
    
    //if(DiffHist){
    
    TH1F *DifDoTr2ANdaxvAzim = new TH1F("DifDoTr2AndaxvAzim","Diference between Azimutal angles",360,-180,180);
    TH1F *DifDoTr2ANdaxvZenit = new TH1F("DifDoTr2AndaxvZenit","Diference between Zenit angles",360,-180,180);
    
    //}
    
    /*************************************************************************************/
    //Histogramm for Pamela's angular efficiency
    /*************************************************************************************/
    
    TH2F *PamAngEffhigh = new TH2F("PamAngEffhigh","Zenit Angle, deg; Azimutal Angle, deg",360,0,360,100,80,180);
    TH2F *PamAngEfflow = new TH2F("PamAngEfflow","Zenit Angle, deg; Azimutal Angle, deg",360,0,360,100,80,180);
    TCanvas *c7 = new TCanvas("c7","Pamela's angular efficiency", 1600,2000);
    
    /*************************************************************************************/
    //Hystogramm to count how much time Pamela have each PitchAngle
    /*************************************************************************************/
    
    TH1F *PitchTime = new TH1F("PitchTime","Time for Pitch",360,0,360);
    TCanvas *c8 = new TCanvas("c8","PitchTime",1200,600);
    
    TH1F *PitchExposure = new TH1F("PitchExposure","Time for Pitch",360,0,360);
    TCanvas *c9 = new TCanvas("c9","PitchTime",1200,600);
    
    /***************************************************************************************/
    //Initialization histogramms for count of protons depends from Energy. Volodia.
    /***************************************************************************************/
/*    
    TH1F* proton10log= new TH1F("Protons Flux","",80,-1.,3.);
    proton10log->GetXaxis()->SetTitle("log10(Energy) , GeV");
    TH1F* proton1log= new TH1F("Protons Flux","",80,-1.,3.);
    TH1F* binw= new TH1F("w","",80,-1.,3.);
    proton1log->GetXaxis()->SetTitle("log10(Energy) , GeV");
    proton1log->SetLineColor(kRed);
    TCanvas *plogCanvas = new TCanvas("proton flux","protonflux", 800, 800);
    // this is bins wide to calculate Flux
    Float_t x;
    Float_t xmin, xmax;
    Float_t binwide[80];
    for (Int_t l=0 ; l<80; l++) {
	binw->Fill(-1.+0.05*l,pow(10.,0.05*(Float_t)(l+1)-1.)- pow(10.,0.05*(Float_t)(l)-1.)); 
    }
   
    for (Int_t l=0 ; l<80; l++) {
    binwide[l]=binw->GetBinContent(l+1);
    }
*/    
    /****************************************************************************/
    //Geometrical Factor. Volodia
    /****************************************************************************/
/*
    TH1F* Geomfactor = new TH1F("G","G",1000,0.1,500.1);
    TH1F* Geomfactorlog= new TH1F("Glog","Glog",80,-1.,3.);
    TH1F* Geomfactorlogelectron= new TH1F("Gloge","Gloge",80,-1.,3.);
    for (Int_t l=0 ; l <80 ;l++) {
	x=pow(10.,(Float_t) 0.05*l-1.);
	Geomfactorlogelectron->SetBinContent(l+1,pow(1/pow(66.7051+50.05404*log10(x),3.5083)+1./pow(21.6238,3.5083),(-1./3.5083)));
	x=ConvertT2R(x,0.938,1., 1.);
	Geomfactorlog->SetBinContent(l+1,pow(1/pow(66.7051+50.05404*log10(x),3.5083)+1./pow(21.6238,3.5083),(-1./3.5083)));
    }//geomfactor for linear scale 
    for (Int_t l=0 ; l <1000 ;l++) {
	x=(Float_t) 0.5*l+0.1;
	Geomfactor->SetBinContent(l+1,pow(1/pow(66.7051+50.05404*log10(x),3.5083)+1./pow(21.6238,3.5083),(-1./3.5083)));
    }
*/
    /****************************************************************************/
    //General loop
    /****************************************************************************/

    for(ULong_t i=0; i<nevents; i++){
    
	T->GetEntry(i);
	
	if(pam_events->GetTrkLevel2()->GetNTracks()==1){
	    
	    /*************************************************************************/
	    //Getting general parameters
	    /*************************************************************************/
	    
	    track = pam_events->GetTrack(0);
	    Int_t M0DE = pam_events->GetOrbitalInfo()->mode;  //0 is zero here
	    Float_t Bx = -pam_events->GetOrbitalInfo()->Bdown;
	    Float_t By = pam_events->GetOrbitalInfo()->Beast;
	    Float_t Bz = pam_events->GetOrbitalInfo()->Bnorth;
	    Float_t Babs = pam_events->GetOrbitalInfo()->Babs;
	    Float_t L = pam_events->GetOrbitalInfo()->L;
	    Float_t rigev = 1/track->GetTrkTrack()->GetDeflection();
	    
	    /*************************************************************************/
	    //Track selection
	    /*************************************************************************/
	    
	    ntr = 0;
	    if ((track->GetTrkTrack()->chi2<=0) || (track->GetTrkTrack()->chi2>=500)) ntr=-1;
	    for(Int_t ii=1;ii<=12;ii++){
		if ((track->GetToFTrack()->beta[ii]<0) || (track->GetToFTrack()->beta[ii]>99)) ntr=-1;
	    }
	    if((track->GetTrkTrack()->GetNX()<=4)&&(track->GetTrkTrack()->GetNY()<4)) ntr=-1;
	    if((M0DE!=0)&&(M0DE!=1)&&(M0DE!=6)&&(M0DE!=2)&&(M0DE!=3)&&(M0DE!=8)&&(M0DE!=4)) ntr=-1;
	    if(sqrt(pow(Bx,2)+pow(By,2)+pow(Bz,2))<0.001) ntr=-1;
	    //if(rigev<0.600) ntr=-1;
	    
	    if(ntr!=-1 || PamExp){
	    
		/*******************************************************************************/
		//Transit reference system
	        /*******************************************************************************/
	    
		TMatrixD Fij = PO->ECItoGreenwich(PO->QuatoECI(pam_events->GetOrbitalInfo()->q0,pam_events->GetOrbitalInfo()->q1,pam_events->GetOrbitalInfo()->q2,pam_events->GetOrbitalInfo()->q3),pam_events->GetOrbitalInfo()->absTime);
		TMatrixD Dij = PO->GreenwichtoGEO(pam_events->GetOrbitalInfo()->lat,pam_events->GetOrbitalInfo()->lon,Fij);
		TMatrixD Iij = PO->ColPermutation(Dij);
	    
	        /*****************************************************************************/
		//Calculate Zenit and Azimutal angle and vector of particle flight in Pamela
	        //using axv[0] and ayv[0] variables
		/*****************************************************************************/	    
	    
		if(!DoTr || DiffHist){
	    
	        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);
	    
		
		}
	    
	        /*****************************************************************************/
		//Calculate Zenit and Azimutal angle and vector of particle flight in Pamela
	        //using DoTrack2 procedure
		/*****************************************************************************/
	    
		if(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.
		/******************************************************************************/
/*		
		if(14.9/L/L > 0. && 14.9/L/L < 1. & Babs >0.24) {
    		Int_t detr = 0;
    		//if((rigev < -0.8  && betaev >0.8)||(rigev>-0.8 && rigev < 0.&& betaev > 0.8) ) //electron2->Fill(log10(fabs(rigev)));
		    if(detr  < (15*exp(-1.6*(fabs(rigev)-0.5))+4)  && rigev >0.5  || (detr < 30 && fabs(rigev) <0.5)) {
			Float_t Ekin=ConvertR2T(rigev,0.938,1., 1.);
			proton1log->Fill(log10(Ekin));
		    }
		}
*/		
		/************/
		//Exposure
		/************/
		
		//for (Int_t itime=0; itime<16;itime++){
		//    if(14.9/L/L > itime && 14.9/L/L < (itime+1) & Babs >0.24) exposure[itime]=exposure[itime]+livetime/1000.;}
		
		/***********************************************************************/
		//Polar area. To calculate angular efficiencÐÐy of registration
		/***********************************************************************/
		
		if(PamExp){

		Float_t AnglTime = 0.16*pam_events->GetTrigLevel2()->dltime[0]+0.01*pam_events->GetTrigLevel2()->dltime[1];
		for(Int_t i=0;i<=360;i++){
		    for(Int_t j=0;j<=100;j++){
		        if(hhigh->GetBin(i,j)>0){
		    	    Double_t pamZenit=(j+80-0.5)*TMath::DegToRad();
		    	    Double_t pamAzim=(i*-0.5)*TMath::DegToRad();
		    	    Double_t _z = cos(pamZenit);
			    Double_t _x = sin(pamZenit)*cos(pamAzim);
		    	    Double_t _y = sin(pamZenit)*sin(pamAzim);
			    TMatrixD Uij = PO->PamelatoGEO(Iij,_x,_y,_z);
		    	    Double_t SP = PO->GetPitchAngle(Uij(0,0),Uij(1,0),Uij(2,0),Bx,By,Bz);
			    PitchExposure->Fill(SP,AnglTime);
 			}
		    }
		}
		
		}
		
		if(PamEff){
		
		if(L>6){
		Float_t Mass = 0;
    		Int_t detr = 0;
    		Float_t Ekin = 0;
    		Float_t betaev = 0.25*(track->GetToFTrack()->beta[0]+track->GetToFTrack()->beta[1]+track->GetToFTrack()->beta[2]+track->GetToFTrack()->beta[3]);    
		    if(14.9/L/L > 0. && 14.9/L/L < 1. & Babs >0.24) {	
    			if((rigev < -0.8  && betaev >0.8)||(rigev>-0.8 && rigev < 0.&& betaev > 0.8) ) {
    			    Ekin=ConvertR2T(rigev,0.511,-1.);
    			    if (Ekin>1) PamAngEffhigh->Fill(TMath::RadToDeg()*PamAzim,TMath::RadToDeg()*PamZenit);
    			    if (Ekin<1) PamAngEfflow->Fill(TMath::RadToDeg()*PamAzim,TMath::RadToDeg()*PamZenit);
    			}//electron2->Fill(log10(fabs(rigev)));
			
			if(detr  < (15*exp(-1.6*(fabs(rigev)-0.5))+4)  && rigev >0.5  || (detr < 30 && fabs(rigev) <0.5)) {
			    Ekin=ConvertR2T(rigev,0.938,-1.);
			    if (Ekin>1) PamAngEffhigh->Fill(TMath::RadToDeg()*PamAzim,TMath::RadToDeg()*PamZenit);
			    if (Ekin<1) PamAngEfflow->Fill(TMath::RadToDeg()*PamAzim,TMath::RadToDeg()*PamZenit);
			}
		    }
		}
		
		}
		
		/***********************************************************************/
		//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){

		/***********************************************************************/
		//Filling angular histogramms for Brazil
		/***********************************************************************/

		if((Babs<0.19)&&(L<1.75)){
			Argv = 0.16*pam_events->GetTrigLevel2()->dltime[0];
			PitchExpositionBrazil->Fill(SB,Argv);
			MainAxisPamelaExpositionBrazil->Fill(SA,Argv);
			PitchExpositionNoOrientationBrazil->Fill(SC,Argv);
			PitchCountBrazil->Fill(SB);
			MainAxisPamelaCountBrazil->Fill(SA);
			PitchCountNoOrientationBrazil->Fill(SC);
			ZenitAngle->Fill(ZC);
			PitchvsLatInBrazil->Fill(pam_events->GetOrbitalInfo()->lat,SA);
			LieveTimevsLatInBrazil->Fill(pam_events->GetOrbitalInfo()->lat,Argv);
		}
		
		/***********************************************************************/
		//Filling angulars histogramms for areas without Brazil
		/***********************************************************************/
		
		if(Babs>0.24){
			Argv = 0.16*pam_events->GetTrigLevel2()->dltime[0];
			PitchExpositionWithoutBrazil->Fill(SB,Argv);
			MainAxisPamelaExpositionWithoutBrazil->Fill(SA,Argv);
			PitchExpositionNoOrientationWithoutBrazil->Fill(SC,Argv);
			PitchCountWithoutBrazil->Fill(SB);
			MainAxisPamelaCountWithoutBrazil->Fill(SA);
			PitchCountNoOrientationWithoutBrazil->Fill(SC);
			ZenitAngle->Fill(ZC);
			
			/***************************************************************/
			//Filling angulars histogramm for equator
			/***************************************************************/
			
			if(L<1.2){
			    Argv = 0.16*pam_events->GetTrigLevel2()->dltime[0];
			    PitchExpositionEquator->Fill(SB,Argv);
			    MainAxisPamelaExpositionEquator->Fill(SA,Argv);
			    PitchExpositionNoOrientationEquator->Fill(SC,Argv);
			    PitchCountEquator->Fill(SB);
			    MainAxisPamelaCountEquator->Fill(SA);
			    PitchCountNoOrientationEquator->Fill(SC);
			    ZenitAngle->Fill(ZC);
			}
		}
		
		} // if(AnglHist)
	    
	    } // if(PamEff||AnglHist)
	    
	    } // if ntr!=-1;
	    
	} // if GetNTrack==1;
    
    } //general loop

    Double_t BinContent;
    
    /*********************************************************/
    //Divide canvases for angular histogramms
    /*********************************************************/
    
    if(AnglHist){
    
    c1->Divide(1,6);
    c2->Divide(1,6);
    c3->Divide(1,6);
    c4->Divide(1,7);
    c5->Divide(1,2);
    
    }
    
    /***************************************************************************************/
    //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){
    
    c7->Divide(1,2);
    
    }
    
    /**************************************************************************************/
    //Draw Angular Histogramms
    /**************************************************************************************/

    if(AnglHist){

    c1->cd(1);
    MainAxisPamelaExpositionBrazil->Draw();
    c1->cd(2);
    MainAxisPamelaCountBrazil->Draw();
    c1->cd(3);
    PitchExpositionNoOrientationBrazil->Draw();
    c1->cd(4);
    PitchCountNoOrientationBrazil->Draw();
    c1->cd(5);
    DiferenceBetweenNoOrientationAndWithOrientationBrazil->Divide(PitchExpositionNoOrientationBrazil,MainAxisPamelaExpositionBrazil);
    DiferenceBetweenNoOrientationAndWithOrientationBrazil->Draw();
    c1->cd(6);
    for(Int_t iu=0;iu<360;iu++){
	BinContent = DiferenceBetweenNoOrientationAndWithOrientationBrazil->GetBinContent(iu);
	DifferenceCountInBrazil->Fill(BinContent);
    }
    DifferenceCountInBrazil->Draw();
    //DifferenceCountInBrazil->SaveAs((TString)outdir+"PIC001A_BrazilExpositionAndCount.root");
    c1->SaveAs((TString)outdir+"PIC001A_BrazilExpositionAndCount.jpg");
    
    c2->cd(1);
    MainAxisPamelaExpositionWithoutBrazil->Draw();
    c2->cd(2);
    MainAxisPamelaCountWithoutBrazil->Draw();
    c2->cd(3);
    PitchExpositionNoOrientationWithoutBrazil->Draw();
    c2->cd(4);
    PitchCountNoOrientationWithoutBrazil->Draw();
    c2->cd(5);
    DiferenceBetweenNoOrientationAndWithOrientationWithoutBrazil->Divide(PitchExpositionNoOrientationWithoutBrazil,MainAxisPamelaExpositionWithoutBrazil);
    DiferenceBetweenNoOrientationAndWithOrientationWithoutBrazil->Draw();
    c2->cd(6);
    for(Int_t iu=0;iu<360;iu++){
	BinContent = DiferenceBetweenNoOrientationAndWithOrientationWithoutBrazil->GetBinContent(iu);
	DifferenceCountWithoutBrazil->Fill(BinContent);
    }
    DifferenceCountWithoutBrazil->Draw();
    c2->SaveAs((TString)outdir+"PIC003A_WithoutBrazilExpositionAndCount.jpg");
    
    c3->cd(1);
    MainAxisPamelaExpositionEquator->Draw();
    c3->cd(2);
    MainAxisPamelaCountEquator->Draw();
    c3->cd(3);
    PitchExpositionNoOrientationEquator->Draw();
    c3->cd(4);
    PitchCountNoOrientationEquator->Draw();
    c3->cd(5);
    DiferenceBetweenNoOrientationAndWithOrientationEquator->Divide(PitchExpositionNoOrientationEquator,MainAxisPamelaExpositionEquator);
    DiferenceBetweenNoOrientationAndWithOrientationEquator->Draw();
    c3->cd(6);
    for(Int_t iu=0;iu<360;iu++){
	BinContent = DiferenceBetweenNoOrientationAndWithOrientationEquator->GetBinContent(iu);
	DifferenceCountEquator->Fill(BinContent);
    }
    DifferenceCountEquator->Draw();
    c3->SaveAs((TString)outdir+"PIC005A_EquatorExpositionAndCount.jpg");
    
    c4->cd(1);
    PitchExpositionBrazil->Draw();
    c4->cd(2);
    PitchCountBrazil->Draw();
    c4->cd(3);
    PitchExpositionWithoutBrazil->Draw();
    c4->cd(4);
    PitchCountWithoutBrazil->Draw();
    c4->cd(5);
    PitchExpositionEquator->Draw();
    c4->cd(6);
    PitchCountEquator->Draw();
    c4->cd(7);
    ZenitAngle->Draw();
    c4->SaveAs((TString)outdir+"PIC006A_PitchCountsAndExpositions.jpg");

    c5->cd(1);
    PitchvsLatInBrazil->Draw();
    c5->cd(2);
    LieveTimevsLatInBrazil->Draw();
    c5->SaveAs((TString)outdir+"PIC006A_PitchvsLotitudeBrazil.jpg");

    }

    /***************************************************************************************/
    //Draw hystogramms for diferences in results between calculating using DoTrack2
    //and using axv[0], ayv[0]
    /**************************************************************************************/

    if(DiffHist){

    c6->cd(1);
    DifDoTr2ANdaxvAzim->Draw();
    c6->cd(2);
    DifDoTr2ANdaxvZenit->Draw();
    c6->SaveAs((TString)outdir+"PIC001Diference.jpg");
    
    }
    
    /***********************************************************************/
    //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
    /**************************************************************************************/

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

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

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

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

    proton1log->Draw("");

*/
}