/[PAMELA software]/ParDirCalc/QtoInclination.C
ViewVC logotype

Diff of /ParDirCalc/QtoInclination.C

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.2 by pam-mep, Mon Jan 21 13:06:43 2008 UTC revision 1.3 by pam-mep, Tue Sep 16 09:38:49 2008 UTC
# Line 6  Line 6 
6  //  //
7  #include "TrAng.h"  #include "TrAng.h"
8  #include "QtoInc.h"  #include "QtoInc.h"
9    #include "option.h"
10    
11  #include <fstream>  #include <fstream>
12  #include <iostream>  #include <iostream>
# Line 24  Line 25 
25  #include <TCanvas.h>  #include <TCanvas.h>
26  //#include <TMatrixD.h>  //#include <TMatrixD.h>
27  #include "TStyle.h"  #include "TStyle.h"
28  #include <TH2F.h>  //#include <TH2F.h>
29    #include <TGraphErrors.h>
30    #include <TMultiGraph.h>
31    
32  #include <PamLevel2.h>  #include <PamLevel2.h>
33  #include <TStyle.h>  #include <TStyle.h>
34    
35  using namespace std;  using namespace std;
36    
37    TH2F* NormalizeTH2F(TH2F* h1,Int_t Nev){
38        TH2F* h2 = new TH2F("h2","Zenit, deg; Azimut, deg",36,0,360,20,140,180);
39        //h2-SetName("h2");
40        //h2.SetBins(36,0,360,20,140,180);
41        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);
42        return h2;
43    }
44    
45    void DivideHystTH2F(TString file1, TString file2, TString hyst1, TString hyst2, Int_t Nev1, Int_t Nev2, string outdir,Bool_t Err){
46        gStyle->SetPalette(1);
47        TFile* f1 = new TFile(file1);
48        TFile* f2 = new TFile(file2);
49        TH2F* h1 = (TH2F*)f1->Get(hyst1);
50        TH2F* h2 = (TH2F*)f2->Get(hyst2);
51        TH2F* h1n = NormalizeTH2F(h1,1);
52        TH2F* h2n = NormalizeTH2F(h2,1);
53        /*TH2F* h1n = new TH2F("h1n","Zenit, deg; Azimut, deg",36,0,360,20,140,180);
54        TH2F* h2n = new TH2F("h2n","Zenit, deg; Azimut, deg",36,0,360,20,140,180);
55        for(Int_t i=0;i<=36;i++) for(Int_t j=0;j<=20;j++){
56            h1n->Fill(i*10-5,j*2+140,h1->GetBinContent(i,j)/h1->GetEntries());
57            h2n->Fill(i*10-5,j*2+140,h2->GetBinContent(i,j)/h2->GetEntries());
58        }*/
59        //vector<TGraphErrors> GrafErr(8);
60        TH2F* h3 = new TH2F("h3","Zenit angle, deg; Azimut angle, deg",36,0,360,20,140,180);
61        h3->Divide(h1n,h2n);
62        TH2F* h3ErrPer;
63        if(Err){
64            TH2F* h3Err = new TH2F("h3Err","Errot of Zenit angle, deg; Azimut angle, deg",36,0,360,20,140,180);
65            h3ErrPer = new TH2F("h3Err","Errot of Zenit angle, deg; Azimut angle, deg",36,0,360,20,140,180);
66            //ofstream fe((TString)outdir+"Divide_errors"+hyst1+"_"+hyst2+".txt");
67            ofstream fe2((TString)outdir+"Protons_Divide_errors"+hyst1+"_"+hyst2+".txt");
68            //Int_t Nev1 = 0;
69            //Int_t Nev2 = 0;
70            //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);}
71            //cout<<"Nev1 = "<<Nev1<<" Nev2 = "<<Nev2<<endl;
72            
73            for(Int_t j=20;j>=1;j--){
74                for(Int_t i=1;i<=36;i++){
75                    //fe<<i<<"\t"<<j<<"\t"<<h1->GetBinContent(i,j)<<"/"<<h2->GetBinContent(i,j)<<" = ";
76                    if(h2->GetBinContent(i,j)!=0){
77                        Float_t Val = (h1->GetBinContent(i,j))/(h2->GetBinContent(i,j));
78                        Float_t Err = 0;
79                        if(h1->GetBinContent(i,j)==0)Err=0;else Err = Val*(sqrt(1./h1->GetBinContent(i,j)+1./h2->GetBinContent(i,j)));
80                        //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);
81                        //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);
82                        Float_t ErrPer = Err*100/(fabs(Val-1));
83                        fe2<<(-j+20)*2+1<<"\t"<<i*10-5<<"\t"<<h1->GetBinContent(i,j)<<"\t"<<h2->GetBinContent(i,j)<<endl;
84                        //fe<<Val<<" +/- "<<Err<<"\t";
85                        if(ErrPer>100)ErrPer=100;
86                        h3ErrPer->Fill(i*10-5,j*2+139,ErrPer);
87                        h3Err->Fill(i*10-5,j*2+139,Err);
88                    }
89                }
90            //fe<<"\n";
91            }
92    
93            //h1ErrCanv->Divide(1,8);
94            for(Int_t i=0;i<8;i++){
95                Float_t x[36];
96                Float_t y[36];
97                Float_t y1[36];
98                Float_t Ex[36];
99                Float_t Ey[36];
100                TH1F* Profil1D = new TH1F("Profil1D","Profil1D",36,0,36);
101                Profil1D->SetFillColor(2);
102                for(Int_t j=0;j<36;j++){
103                    y[j] = h3->GetBinContent(j,19-i);
104                    Ey[j] = h3Err->GetBinContent(j,19-i);
105                    x[j] = j;
106                    Ex[j] = 0;
107                    y1[j] = 1;
108                    Profil1D->Fill(j,h2->GetBinContent(j,19-i));
109                }
110                TCanvas* h1ErrCanv = new TCanvas("h1ErrCanv","",1200,800);
111                h1ErrCanv->Divide(1,2);
112                TGraphErrors* Gep = new TGraphErrors(36,x,y,Ex,Ey);
113                TGraph* _1 = new TGraph(36,x,y1);
114                Gep->SetLineColor(4);
115                Gep->SetFillColor(2);
116                h1ErrCanv->cd(1);
117                Gep->Draw("AB");
118                _1->Draw("");
119                h1ErrCanv->cd(2);
120                Profil1D->Draw("");
121                /*string s;
122                stringstream out;
123                out<<i;
124                s=out.str();
125                h1ErrCanv->SaveAs((TString)outdir+"Divide_"+hyst1+"_"+hyst2+"_1Ds"+s+".png");*/
126                Gep->Delete();
127                delete h1ErrCanv;
128            }
129            /*fe2.close();
130            TFile f3((TString)outdir+"DivideHyst"+hyst1+"_"+hyst2+".root","recreate");
131            h3->Write();
132            f3.Close();*/
133        }    
134        /*TCanvas* c1 = new TCanvas("c1",hyst1+"/"+hyst2,1200,1200);
135        if(Err) c1->Divide(1,2);
136        c1->cd(1);
137        h3->Draw("colz");
138        if(Err){
139            c1->cd(2);
140            h3ErrPer->Draw("colz");
141        }
142        c1->SaveAs((TString)outdir+"Divide_"+hyst1+"_"+hyst2+".png");
143        delete c1;*/
144        h2->Delete();
145        h1->Delete();
146        h3->Delete();
147        f1->Delete();
148        f2->Delete();
149    }
150    
151    void DiferenceHystTH2F(TString file1, TString file2, TString hyst1, TString hyst2, Int_t Nev1, Int_t Nev2, string outdir){
152        gStyle->SetPalette(1);
153        TFile* f1 = new TFile(file1);
154        TFile* f2 = new TFile(file2);
155        TH2F* h1 = (TH2F*)f1->Get(hyst1);
156        TH2F* h2 = (TH2F*)f2->Get(hyst2);
157        TH2F* h1n = NormalizeTH2F(h1,Nev1);
158        TH2F* h2n = NormalizeTH2F(h2,Nev2);
159        TH2F* h3 = new TH2F("h3","Zenit angle, deg; Azimut angle, deg",36,0,360,20,140,180);
160        TH2F* h3Err = new TH2F("h3Err","Zenit angle, deg; Azimut angle, deg",36,0,360,20,140,180);
161        TH2F* h3ErrPer = new TH2F("h3Err","Errot of Zenit angle, deg; Azimut angle, deg",36,0,360,20,140,180);
162        ofstream fe((TString)outdir+"Diference_errors2.txt");
163        for(Int_t i=1;i<=36;i++){
164            for(Int_t j=1;j<=20;j++){
165                Float_t Val = h1->GetBinContent(i,j)/Nev1-h2->GetBinContent(i,j)/Nev2;
166                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));
167                Float_t ErrPer = Err*100/(fabs(Val));
168                h3->Fill(i*10-5,j*2+139,h1->GetBinContent(i,j)/Nev1-h2->GetBinContent(i,j)/Nev2);
169                fe<<i<<"\t"<<j<<"\t"<<h1->GetBinContent(i,j)/Nev1<<" - "<<h2->GetBinContent(i,j)/Nev2<<" = "<<Val;
170                if(h1->GetBinContent(i,j)!=0 && h2->GetBinContent(i,j)!=0){
171                    fe<<" +/- "<<Err<<" (+/- "<<ErrPer<<" )\n";
172                    if(ErrPer>100)ErrPer=100;
173                    h3Err->Fill(i*10-5,j*2+139,Err);
174                    h3ErrPer->Fill(i*10-5,j*2+139,ErrPer);
175                }else fe<<"\n";
176            }
177        fe<<"\n";
178        }
179        for(Int_t i=0;i<8;i++){
180            Float_t x[36];
181            Float_t y[36];
182            Float_t Ex[36];
183            Float_t Ey[36];
184            TH1F* Profil1D = new TH1F("Profil1D","Profil1D",36,0,36);
185            Profil1D->SetFillColor(2);
186            for(Int_t j=0;j<36;j++){
187                if(h1n->GetBinContent(j,19-i)!=0)y[j] = h3->GetBinContent(j,19-i)*100/h1n->GetBinContent(j,19-i); else y[j]=0;
188                if(h1n->GetBinContent(j,19-i)!=0) Ey[j] = h3Err->GetBinContent(j,19-i)*100/h1n->GetBinContent(j,19-i); else Ey[j]=0;
189                Profil1D->Fill(j,h2->GetBinContent(j,19-i));
190                //cout<<"h3 = "<<h3->GetBinContent(j,19-i)<<" h1n = "<<h1n->GetBinContent(j,19-i)<<" yn = "<<y[j]<<endl;
191                //Ey[j] = h3ErrPer->GetBinContent(j,19-i);
192                //cout<<Ey[j]<<endl;
193                x[j] = j;
194                Ex[j] = 0;
195            }
196            TCanvas* h1ErrCanv = new TCanvas("h1ErrCanv","",1200,800);
197            h1ErrCanv->Divide(1,2);
198            TGraphErrors* Gep = new TGraphErrors(36,x,y,Ex,Ey);
199            Gep->SetLineColor(4);
200            Gep->SetFillColor(2);
201            h1ErrCanv->cd(1);
202            Gep->Draw("AB");
203            h1ErrCanv->cd(2);
204            Profil1D->Draw("");
205            string s;
206            stringstream out;
207            out<<i;
208            s=out.str();
209            h1ErrCanv->SaveAs((TString)outdir+"Diference_"+hyst1+"_"+hyst2+"_1Ds"+s+".png");
210            Gep->Delete();
211            Profil1D->Delete();
212            delete h1ErrCanv;
213        }
214        fe.close();
215        TFile f3((TString)outdir+"resultDifferenceHyst.root","recreate");
216        h3->Write();
217        f3.Close();
218        TCanvas* c1 = new TCanvas("c1",hyst1+" - "+hyst2,1200,1200);
219        c1->Divide(1,2);
220        c1->cd(1);
221        h3->Draw("colz");
222        c1->cd(2);
223        h3ErrPer->Draw("colz");
224        c1->SaveAs((TString)outdir+"Difference_"+hyst1+"_"+hyst2+".png");
225        delete c1;
226        h2->Delete();
227        h1->Delete();
228        h3->Delete();
229        f1->Delete();
230        f2->Delete();
231    }
232    
233    void MergeHystTH2F(TString file1, TString file2, TString hyst1, TString hyst2, string outdir){
234        gStyle->SetPalette(1);
235        TFile* f1 = new TFile(file1);
236        TFile* f2 = new TFile(file2);
237        TH2F* h1 = (TH2F*)f1->Get(hyst1);
238        TH2F* h2 = (TH2F*)f2->Get(hyst2);
239        TH2F* h3 = new TH2F("h3","Zenit angle, deg; Azimut angle, deg",36,0,360,20,140,180);
240        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));
241        TFile f3((TString)outdir+hyst1+"_"+hyst2+"_resultMergeHyst.root","recreate");
242        h3->Write();
243        f3.Close();
244        TCanvas* c1 = new TCanvas("c1",hyst1+" + "+hyst2,1200,800);
245        c1->cd();
246        h3->Draw("colz");
247        c1->SaveAs((TString)outdir+"Merge_"+hyst1+"_"+hyst2+".png");
248        delete c1;
249        h2->Delete();
250        h1->Delete();
251        h3->Delete();
252        f1->Delete();
253        f2->Delete();
254    }
255    
256    
257    void Convert2Dto1Dhyst(TString file, TString hyst, string outdir){
258        TFile* f1 = new TFile(file);
259        TH2F* h1 = (TH2F*)f1->Get(hyst);
260        //TH1F* h2 = new TH1F("h2","Zenit Efficiency",20,140,180);
261        //TH1F* h3 = new TH1F("h3","Azimuth Efficiency",36,0,360);
262        //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));
263        //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));
264        TH1D* h2 = (TH1D*)h1->ProjectionY(" ",1,20,"");
265        TCanvas* c1 = new TCanvas("c1","Zenit Efficiency",1200,800);
266        c1->cd();
267        h2->Draw();
268        c1->SaveAs((TString)outdir+hyst+"Zenit.png");
269        delete c1;
270        TFile f2((TString)outdir+hyst+"Zenit.root","recreate");
271        h2->Write();
272        f2.Close();
273        TH1D* h3 = (TH1D*)h1->ProjectionX(" ",1,36,"");
274        TCanvas* c2 = new TCanvas("c2","Azimuth Efficiency",1200,800);
275        c2->cd();
276        h3->Draw();
277        c2->SaveAs((TString)outdir+hyst+"Azimuth.png");
278        delete c2;
279        TFile f3((TString)outdir+hyst+"Azimuth.root","recreate");
280        h3->Write();
281        f3.Close();
282        h3->Delete();
283        h1->Delete();
284        h2->Delete();
285        f1->Delete();
286    }
287    
288    void ReBinTH2F(TString file, TString hyst, string outdir){
289        gStyle->SetPalette(1);
290        TFile* f1 = new TFile(file);
291        TH2F* h1 = (TH2F*)f1->Get(hyst);
292        TH2F* h2 = new TH2F("h2","Zenith, deg; Azimuth, deg",18,0,360,20,140,180);
293        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++;}
294        TFile f2((TString)outdir+hyst+"Re-Bin.root","recreate");
295        h2->Write();
296        f2.Close();
297        TCanvas* c1 = new TCanvas("c1",hyst,800,600);
298        c1->cd();
299        h2->Draw("colz");
300        c1->SaveAs((TString)outdir+hyst+"Re-Bin.png");
301        delete c1;
302        h2->Delete();
303        h1->Delete();
304        f1->Delete();
305    }
306    
307    void Modeling(Double_t R, Int_t Nevents, Float_t zmax, string outdir){
308        cout<<"Modelling..."<<endl;
309        gStyle->SetPalette(1);
310        TH2F* h1 = new TH2F("h1","Azimut,deg;Zenit,deg",36,0,360,20,140,180);
311        TH2F* h2 = new TH2F("h2","Azimut,deg;Zenit,deg",36,0,360,20,140,180);
312        TH2F* h3 = new TH2F("h3","Azimut,deg;Zenit,deg",36,0,360,20,140,180);
313        TH2F* h4 = new TH2F("h4","Azimut,deg;Zenit,deg",36,0,360,20,140,180);
314        Float_t A;
315        Float_t Zenit;
316        Float_t Azimuth;
317        for(Int_t i=0;i<Nevents;i++){
318            A=rand() % 1000; A=A/1000;
319            Zenit=TMath::RadToDeg()*acos(pow(A,0.5));
320            //Zenit=3;
321            Azimuth=rand() % 360000;
322            Azimuth=Azimuth/1000;
323            h4->Fill(Azimuth,180-Zenit);
324            /*for(Int_t nx=0;nx<24;nx++){
325                for(Int_t ny=0;ny<24;ny++){
326                    h2->Fill(Azimuth,180-Zenit);
327                    Float_t x = (nx+0.5)*6.75;
328                    Float_t y = (ny+0.5)*5.5;
329                    Float_t k = sin(TMath::DegToRad()*Zenit)*cos(TMath::DegToRad()*Azimuth);
330                    Float_t l = sin(TMath::DegToRad()*Zenit)*sin(TMath::DegToRad()*Azimuth);
331                    Float_t m = cos(TMath::DegToRad()*Zenit);
332                    Float_t xl = x+k*zmax/m;
333                    Float_t yl = y+l*zmax/m;
334                    Float_t hl = sqrt(pow(x-xl,2)+zmax);
335                    Float_t za=0; Float_t xa=0;
336                    Float_t Rk=1000*0.3*R/0.4;
337                    //cout<<"Zenith = "<<Zenit<<" Azimuth = "<<Azimuth<<endl;
338                    //cout<<"x = "<<x<<" y = "<<y<<" Rk = "<<Rk<<" xl = "<<xl<<" yl = "<<yl<<endl;
339                    if(x-xl==0){
340                     za==0;
341                     xa==x-Rk;
342                    }else{
343                        //Float_t ko = 0;
344                        //Float_t koo = 0;
345                        if(x-xl<0) {za = zmax-Rk*cos(atan(zmax/fabs(x-xl)));}
346                        if(x-xl>0) {za = zmax+Rk*cos(atan(zmax/fabs(x-xl)));}
347                        xa = x-Rk*sin(atan(zmax/fabs(x-xl)));
348                    }
349                    //cout<<"xa = "<<xa<<" za = "<<za<<" atan = "<<TMath::RadToDeg()*atan(445.1/fabs(x-xl))<<endl;
350                    Float_t xll;
351                    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;
352                    if(xll<=162 && xll>=0 && yl<=132 && yl>=0) h1->Fill(Azimuth,180-Zenit);
353                    //cout<<"xll = "<<xll<<endl;
354                    //if(xl<=162 && xl>=0 && yl<=132 && yl>=0)cin>>za;
355                }
356            }*/
357        }/*
358        TCanvas* c1 = new TCanvas("c1","Model distribution",1200,800);
359        h1->Draw("colz");
360        c1->SaveAs((TString)outdir+".png");
361        TFile f2((TString)outdir+".root","recreate");
362        h1->Write();
363        f2.Close();
364        */
365        TCanvas* c3 = new TCanvas("c3","Model Isotropic distribution",1200,800);
366        h4->Draw("colz");
367        c3->SaveAs((TString)outdir+"_isotrop.png");
368        TFile f4((TString)outdir+"_isotrop.root","recreate");
369        h4->Write();
370        f4.Close();
371        /*
372        TFile f3((TString)outdir+"Abs.root","recreate");
373        h2->Write();
374        f3.Close();
375        h3->Divide(h2,h1);
376        TCanvas* c2 = new TCanvas("c2","Model distribution",1200,800);
377        h3->Draw("colz");
378        c2->SaveAs((TString)outdir+"Rate.png");
379        delete c1;
380        delete c2;*/
381        delete c3;
382        h1->Delete();
383        h2->Delete();
384        h3->Delete();
385        h4->Delete();
386    }
387            
388  int main(int argc, char* argv[]){  int main(int argc, char* argv[]){
389    
390    OptionParam opt1;
391    
392  string outdir = "/home/pamelaprod/malakhov/QtoInc/Picture/";  string outdir = "/home/pamelaprod/malakhov/QtoInc/Picture/";
393  Bool_t AnglHist = false;  Bool_t AnglHist = false;
394  Bool_t DiffHist = false;  Bool_t DiffHist = false;
# Line 41  Bool_t DoTr = false; Line 397  Bool_t DoTr = false;
397  Bool_t PamAngTime = false;  Bool_t PamAngTime = false;
398  Bool_t PamExp = false;  Bool_t PamExp = false;
399    
400    /*vector<TH1F> Hyst(1);
401    Hyst[0].SetName("ElectronLowEnergy");
402    Hyst[0].SetBins(36,0,360);
403    Hyst[0].FillRandom("gaus",5000);
404    TFile f1("/home/pamelaprod/malakhov/QtoInc/TEST/TEST.root","recreate");
405    Hyst[0].Write();
406    f1.Close();
407    //TFile* f1 = new TFile("/home/pamelaprod/malakhov/QtoInc/TEST/TEST.root");
408    //TH1F* h1 = (TH1F*)f1->Get("");
409    //cout<<h1->GetEntries()<<endl;
410    */
411    
412  if(argc<2){  if(argc<2){
413      cout<<"You have to insert at least a file to analisys \nUsing --help for more information\n";      cout<<"You have to insert at least a file to analisys \nUsing --help for more information\n";
414      exit(1);      exit(1);
415  }  }
416    //opt1.opt.PE.Proton.Hyst1.resize(opt1.opt.PE.Proton.Hyst1.size()+1);
417    //opt1.opt.PE.Proton.Hyst1[0].SetBins(36,0,360);
418    //opt1.opt.PE.Proton.Hyst1[0].SetName("Test2");
419    //cout<<opt1.opt.PE.Proton.Hyst1.size()<<endl;
420    //cout<<opt1.opt.PE.Proton.Hyst1[0].GetName()<<endl;
421    
422    //vector<TH2F>Hyst2DF(1);
423    //Hyst2DF[0].SetName("Test1");
424    //Hyst2DF[0].SetBins(36,0,360,36,0,360);
425    //Hyst2DF.resize(Hyst2DF.size()+1);
426    //Hyst2DF[1].SetName("Test2");
427    //Hyst2DF[1].SetBins(36,0,360,36,0,360);
428    //cout<<Hyst2DF[0].GetName()<<endl;
429    
430  if(!strcmp(argv[1],"--help")){  if(!strcmp(argv[1],"--help")){
431      cout<<"\nUsage \n";      opt1.helprint();
432    /*    cout<<"\nUsage \n";
433      cout<<"\n QtoInc --WorkDir directory_with_level2_files --OutPuth directory_for_output files [Options] \n";      cout<<"\n QtoInc --WorkDir directory_with_level2_files --OutPuth directory_for_output files [Options] \n";
434      cout<<"\n --WorkDir         full puth to the Level2 file \n";      cout<<"\n --WorkDir         full puth to the Level2 file \n";
435      cout<<" --OutPuth           full puth for putputh files \n";      cout<<" --OutPuth           full puth for putputh files \n";
# Line 61  if(!strcmp(argv[1],"--help")){ Line 443  if(!strcmp(argv[1],"--help")){
443      cout<<"-PamExp              Calculating exposition for each possible direction in Pamela's aperture \n";      cout<<"-PamExp              Calculating exposition for each possible direction in Pamela's aperture \n";
444      cout<<"\nFor Example: \n";      cout<<"\nFor Example: \n";
445      cout<<"\n./QtoInc.exe --WorkDir /data01/_3/704_3/ --OutPuth /home/pamelaprod/malakhov/QtoInc/Picture5/ -DiffHist -DoTr \n\n";      cout<<"\n./QtoInc.exe --WorkDir /data01/_3/704_3/ --OutPuth /home/pamelaprod/malakhov/QtoInc/Picture5/ -DiffHist -DoTr \n\n";
446    */    exit(1);
447    }
448    
449    /*if(!strcmp(argv[1],"--Normalize")){
450        if(4 >= argc){
451            cout<<"--Normalize  needs argument\n See --help\n";
452            exit(1);
453        }
454        NormalizeTH2F((TString)argv[2],(TString)argv[3],argv[4]);
455        exit(0);
456    }*/
457    
458    if(!strcmp(argv[1],"--Model")){
459        if(5 >= argc){
460            cout<<"--Model  needs argument\n See --help\n";
461            exit(1);
462        }
463        Modeling(atof(argv[2]),atoi(argv[3]),atof(argv[4]),argv[5]);
464        exit(0);
465    }
466    
467    if(!strcmp(argv[1],"--Diference")){
468        if(6 >= argc){
469            cout<<"--Diference  needs argument\n See --help\n";
470            exit(1);
471        }
472        
473        DiferenceHystTH2F((TString)argv[2],(TString)argv[3],(TString)argv[4],(TString)argv[5],atoi(argv[6]),atoi(argv[7]),argv[8]);
474        exit(0);
475    }
476    
477    if(!strcmp(argv[1],"--Compare")){
478        //gStyle->SetPalette(1);
479        //TMultiGraph* mg;
480        //vector<TGraphErrors*> Gep;
481        Float_t Val[20][36];
482        Float_t Err[20][36];
483        string tmp;
484        char g;
485        TString out = "";
486        for(Int_t ii=2;ii<argc;ii++){
487            if(!strcmp(argv[ii],"-out")){
488                out = (TString)argv[ii+1];
489                ii++;
490            }else{
491                if(!strcmp(argv[ii],"-in"));
492                ifstream in(argv[ii+1],ios::in);
493                for(Int_t i=20;i>0;i--){
494                for(Int_t j=1;j<=36;j++){
495                    in>>tmp;in>>tmp;in>>g;
496                    if(g=='0'){
497                        Val[i-1][j-1]=0;
498                        Err[i-1][j-1]=0;
499                        Bool_t xx=false;
500                        Bool_t yy=false;
501                        while(g!='='){
502                            in>>g;
503                            if(g=='/'){yy=true;continue;}
504                            if(yy && g=='0') xx=true; else yy=false;
505                        }
506                        if(!xx){
507                            in>>tmp;
508                            in>>tmp;
509                            in>>tmp;
510                        }
511                    }else{
512                        while(g!='/') in>>g;
513                        in>>g;
514                        if(g=='0'){
515                            in>>tmp;
516                            Val[i-1][j-1]=0;
517                            Err[i-1][j-1]=0;
518                        }else{
519                            while(g!='=') in>>g;
520                            in>>Val[i-1][j-1];
521                            in>>tmp;
522                            in>>Err[i-1][j-1];
523                        }
524                    }
525                }}
526                ii++;
527                //Val.resize(Val.size()+1);
528                //Val[Val.size()-1][1][1]=0;
529            }
530        }
531        ofstream fv(out+"Val.txt");
532        ofstream fe(out+"Err.txt");
533        for(Int_t i=20;i>0;i--){
534        for(Int_t j=1;j<=36;j++){fv<<Val[i-1][j-1]<<"\t";fe<<Err[i-1][j-1]<<"\t";}
535            fv<<endl;fe<<endl;
536        }
537        fv.close();
538        fe.close();
539    }
540    
541    if(!strcmp(argv[1],"--Merge")){
542        if(6 >= argc){
543            cout<<"--Merge  needs argument\n See --help\n";
544            exit(1);
545        }
546        gStyle->SetPalette(1);
547        //vector<TString> Hysts;
548        //vector<TString> Files;
549        vector<TH1F*> h1;
550        vector<TH2F*> h2;
551        TH2F h2sum;
552        TH1F h1sum;
553        Bool_t H2 = false;
554        TString out = "";
555        if(!strcmp(argv[2],"2D")) H2=true;
556        for(Int_t i=3;i<argc;i++){
557            if(!strcmp(argv[i],"-out")){
558                out = (TString)argv[i+1];
559                i++;
560                //break;
561            }else{
562                if(H2){
563                    h2.resize(h2.size()+1);
564                    cout<<"h2 resize... New size is "<<h2.size()<<endl;
565                    TFile* f2 = new TFile((TString)argv[i]);
566                    //cout<<"file is reading"<<endl;
567                    //f2->ls();
568                    h2[h2.size()-1] = (TH2F*)f2->Get(argv[i+1]);
569                    //cout<<argv[i+1]<<endl;
570                    //cout<<h3->GetNbinsX()<<endl;
571                    //delete h3;
572                    //delete f1;
573                    i++;
574                }else{
575                    h1.resize(h1.size()+1);
576                    TFile* f1 = new TFile((TString)argv[i]);
577                    h1[h1.size()-1] = (TH1F*)f1->Get(argv[i+1]);
578                    i++;
579                }
580            }
581        }
582        cout<<"outputh is "<<out<<endl;
583        //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;
584        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
585        h1sum.SetBins(h1[0]->GetNbinsX(),h1[0]->GetXaxis()->GetXmin(),h1[0]->GetXaxis()->GetXmax());
586        h1sum.SetName("Proton");
587        h2sum.SetName("Proton");
588        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]);
589        if(H2){
590            TFile fsum(out+".root","recreate");
591            h2sum.Write();
592            fsum.Close();
593            TCanvas* c1 = new TCanvas("c1","resultMerge",1200,800);
594            c1->cd();
595            h2sum.Draw("colz");
596            c1->SaveAs(out+".png");
597            delete c1;
598        }else{
599            TFile fsum(out+".root","recreate");
600            h1sum.Write();
601            fsum.Close();
602            TCanvas* c1 = new TCanvas("c1","resultMerge",1200,800);
603            c1->cd();
604            h1sum.Draw();
605            c1->SaveAs(out+".png");
606            delete c1;
607        }
608        //MergeHystTH2F((TString)argv[2],(TString)argv[3],(TString)argv[4],(TString)argv[5],argv[6]);
609        exit(0);
610    }
611    
612    if(!strcmp(argv[1],"--PlotEfficiency")){
613        vector<TGraphErrors> ge;
614        TString out = "";
615        Int_t t = 1;
616        //TGraphErrors tmp;
617        for(Int_t i=2;i<argc;i++){
618            if(!strcmp(argv[i],"-out")){
619                out = (TString)argv[i+1];
620                i++;
621            }
622            if(!strcmp(argv[i],"-in")){
623                ifstream inVal((TString)argv[i+1]+"Val.txt",ios::in);
624                ifstream inErr((TString)argv[i+1]+"Err.txt",ios::in);
625                Float_t Val[20][36];
626                Float_t Err[20][36];
627                for(Int_t j=0;j<20;j++)for(Int_t k=0;k<36;k++){inVal>>Val[j][k];inErr>>Err[j][k];}
628                if(!strcmp(argv[i+2],"-yline")){
629                    Int_t n = atoi(argv[i+3]);
630                    Float_t x[36];
631                    Float_t y[36];
632                    Float_t xe[36];
633                    Float_t ye[36];
634                    ge.resize(ge.size()+1);
635                    //ge[ge.size()-1]->Set(36);
636                    for(Int_t j=0;j<36;j++){
637                        x[j]=j*10+5;
638                        xe[j]=5;
639                        y[j]=Val[n][j];
640                        ye[j]=Err[n][j];
641                        ge[ge.size()-1].SetPoint(j,x[j],y[j]);
642                        ge[ge.size()-1].SetPointError(j,xe[j],ye[j]);
643                        //tmp.SetPoint(j,x[j],y[j]);
644                    }
645                    ge[ge.size()-1].SetLineColor(t);
646                    ge[ge.size()-1].SetTitle((TString)argv[i+1]+" Pitch  "+(TString)argv[i+3]);
647                    t++;
648                    //TGraphErrors *tmp = new TGraphErrors(36,x,y,xe,ye);
649                    
650                    
651                    //TGraphErrors df;
652                    //df = *tmp;
653                    //cout<<"check"<<endl;
654                    //ge[ge.size()-1] = tmp;
655                    //cout<<i<<endl;
656                    //tmp->Delete();
657                    i+=3;
658                }
659                if(!strcmp(argv[i+2],"-xline")){
660                    Int_t p = 0;
661                    Int_t n=0;
662                    if(!strcmp(argv[i+3],"All")) p=36;else n = atoi(argv[i+3]);
663                    Float_t x[20];
664                    Float_t y[20];
665                    Float_t xe[20];
666                    Float_t ye[20];
667                    for(Int_t j=0;j<20;j++){y[j]=0;ye[j]=0;}
668                    ge.resize(ge.size()+1);
669                    for(Int_t j=0;j<20;j++){
670                        x[j]=j*2+1;
671                        xe[j]=1;
672                        for(Int_t yi=0;yi<p;yi++){
673                            if(p!=0)n=yi;
674                            y[j]+=Val[j][n];
675                            ye[j]+=Err[j][n];
676                        }
677                        ge[ge.size()-1].SetPoint(j,x[j],y[j]);
678                        ge[ge.size()-1].SetPointError(j,xe[j],ye[j]);
679                    }
680                    ge[ge.size()-1].SetLineColor(t);
681                    ge[ge.size()-1].SetTitle((TString)argv[i+1]+" Azimuth  "+(TString)argv[i+3]);
682                    t++;
683                    //TGraphErrors* tmp = new TGraphErrors(20,x,y,xe,ye);
684                    //ge.resize(ge.size()+1);
685                    //ge[ge.size()-1] = (*tmp);
686                    //tmp->Delete();
687                    i+=3;
688                }
689            }
690        }
691        TMultiGraph* mg = new TMultiGraph();
692        cout<<"check"<<endl;
693        for(Int_t i=0;i<ge.size();i++){
694            //ge[i]->SetLineColor(i+1);
695            //cout<<i<<endl;
696            mg->Add(&ge[i]);
697        }
698        TCanvas* c1 = new TCanvas("c1","eryt",1200,800);
699        c1->cd();
700        mg->Draw("AP");
701        cout<<"check"<<endl;
702        cout<<"out = "<<out<<endl;
703        c1->SaveAs(out+".png");
704        delete c1;
705        mg->Delete();
706      exit(1);      exit(1);
707  }  }
708    
709  if(!strcmp(argv[3],"--OutPuth")){  if(!strcmp(argv[1],"--Divide")){
710      if (4 >= argc){      if(argc<7){
711      cout<<"--OutPath needs argument\n See --help\n";              cout<<"--Divide  needs argument\n See --help\n";
712      }else outdir = argv[4];          exit(1);
713        }/*
714        Bool_t H2 = false;
715        TString out = "";
716        if(!strcmp(argv[2],"2D")) H2=true;
717        out = (TString)argv[7];
718        gStyle->SetPalette(1);
719        TFile* f1 = new TFile((TString)argv[3]);
720        TFile* f2 = new TFile((TString)argv[5]);
721        TH1F* h11D; TH1F* h21D;
722        TH2F* h12D; TH2F* h22D;
723        TH1F h1n1D; TH1F h2n1D;
724        TH2F h1n2D; TH2F h2n2D;
725        TH1F h1div; TH2F h2div;
726        cout<<"out is "<<out<<endl;
727        if(H2){
728            h12D = (TH2F*)f1->Get(argv[4]);
729            h22D = (TH2F*)f2->Get(argv[6]);
730            h1n2D.SetBins(h12D->GetNbinsX(),h12D->GetXaxis()->GetXmin(),h12D->GetXaxis()->GetXmax(),h12D->GetNbinsY(),h12D->GetYaxis()->GetXmin(),h12D->GetYaxis()->GetXmax());
731            h2n2D.SetBins(h12D->GetNbinsX(),h12D->GetXaxis()->GetXmin(),h12D->GetXaxis()->GetXmax(),h12D->GetNbinsY(),h12D->GetYaxis()->GetXmin(),h12D->GetYaxis()->GetXmax());
732            h2div.SetBins(h12D->GetNbinsX(),h12D->GetXaxis()->GetXmin(),h12D->GetXaxis()->GetXmax(),h12D->GetNbinsY(),h12D->GetYaxis()->GetXmin(),h12D->GetYaxis()->GetXmax());
733            for(Int_t i = 0; i < h12D->GetNbinsX(); i++){
734                for(Int_t j = 0; j < h12D->GetNbinsY(); j++){
735                    h1n2D.SetBinContent(i,j,h12D->GetBinContent(i,j)/h12D->GetEntries());
736                    h2n2D.SetBinContent(i,j,h22D->GetBinContent(i,j)/h22D->GetEntries());
737                }
738            }
739            h2div.Divide(&h1n2D,&h2n2D);
740        }else{
741            h11D = (TH1F*)f1->Get((TString)argv[4]);
742            h21D = (TH1F*)f2->Get((TString)argv[6]);
743            h1n1D.SetBins(h11D->GetNbinsX(),h11D->GetXaxis()->GetXmin(),h11D->GetXaxis()->GetXmax());
744            h2n1D.SetBins(h21D->GetNbinsX(),h21D->GetXaxis()->GetXmin(),h21D->GetXaxis()->GetXmax());
745            h1div.SetBins(h21D->GetNbinsX(),h21D->GetXaxis()->GetXmin(),h21D->GetXaxis()->GetXmax());
746            for(Int_t i = 0; i<h11D->GetNbinsX();i++){
747                h1n1D.SetBinContent(i,h11D->GetBinContent(i)/h11D->GetEntries());
748                h2n1D.SetBinContent(i,h21D->GetBinContent(i)/h21D->GetEntries());
749            }
750            h1div.Divide(&h1n1D,&h2n1D);
751        }
752        h2div.SetName("Proton");
753        h1div.SetName("Proton");
754        if(H2){
755            TFile fdiv(out+".root","recreate");
756            h2div.Write();
757            fdiv.Close();
758            TCanvas* c1 = new TCanvas("c1","resultDivide",1200,800);
759            c1->cd();
760            h2div.Draw("colz");
761            c1->SaveAs(out+".png");
762            delete c1;
763        }else{
764            TFile fdiv(out+".root","recreate");
765            h1div.Write();
766            fdiv.Close();
767            TCanvas* c1 = new TCanvas("c1","resultDivide",1200,800);
768            c1->cd();
769            h1div.Draw();
770            c1->SaveAs(out+".png");
771            delete c1;
772        }*/
773        Bool_t Err = false;
774        if((argc==10) && !strcmp(argv[9],"-Err")) Err = true;
775        DivideHystTH2F((TString)argv[2],(TString)argv[3],(TString)argv[4],(TString)argv[5],atoi(argv[6]),atoi(argv[7]),argv[8],Err);
776        exit(0);
777  }  }
778    
779  for(Int_t i=4; i<argc; i++){  if(!strcmp(argv[1],"--ReBin")){
780      if(!strcmp(argv[i],"-DoTr")) DoTr = true;      if(4 >= argc){
781      if(!strcmp(argv[i],"-PamEff")) PamEff = true;          cout<<"--ReBin  needs argument\n See --help\n";
782      if(!strcmp(argv[i],"-DiffHist")){DiffHist = true; DoTr = true;}          exit(1);
783      if(!strcmp(argv[i],"-AngHist")) AnglHist = true;      }
784      if(!strcmp(argv[i],"-PamAngTime")) PamAngTime = true;      ReBinTH2F((TString)argv[2],(TString)argv[3],argv[4]);
785      if(!strcmp(argv[i],"-PamExp")) PamExp = true;      exit(0);
786  }  }
787    
788  if(!PamEff && !DiffHist && !AnglHist && !PamAngTime && !PamExp) AnglHist = true;  if(!strcmp(argv[1],"--2Dto1D")){
789        if(4 >= argc){
790            cout<<"--2Dto1D  needs argument\n See --help\n";
791            exit(1);
792        }
793        Convert2Dto1Dhyst((TString)argv[2],(TString)argv[3],argv[4]);
794        exit(0);
795    }
796    
797    if(!strcmp(argv[3],"--OutPuth")){
798        if (4 >= argc){
799        cout<<"--OutPath needs argument\n See --help\n";    
800        }else{
801            outdir = argv[4];
802            opt1.outdir = argv[4];
803        }
804    }else{
805        cout<<"\n ERROR:\tYou need --OutPuth as second argument for run programm\n See --help for more information\n"<<endl;
806        exit(1);
807    }
808    
809  if(!strcmp(argv[1],"--WorkDir")){  if(!strcmp(argv[1],"--WorkDir")){
810      if (2 >= argc){      if (2 >= argc){
811      cout<<"--WorkDir needs argument\n See --help\n";          cout<<"--WorkDir needs argument\n See --help\n";
812            exit(1);
813        }
814    }else{
815        cout<<"\n ERROR:\tYou need --WorkDir as first argument for run programm\n See --help for more information\n"<<endl;
816      exit(1);      exit(1);
     }else Calculate(argv[2],outdir,AnglHist,DiffHist,PamEff,DoTr,PamAngTime,PamExp);  
817  }  }
818    
819    if(argc==5){
820        //opt1.Default();
821        cout<<"MESSAGE:\tSetting default options \n\n";
822    }
823    
824    for(Int_t i=5; i<argc; i++){
825    
826        if(!strcmp(argv[i],"-do")) opt1.debug = true;
827        //opt1.AnglHystFunc(argc,argv);
828        //exit(1);
829        /*
830        Int_t NHyst1D = 0; Int_t NHyst2D = 0;
831        if(!strcmp(argv[i],"--AnglHyst")){
832            
833            if(opt1.debug) cout<<"\n MESSAGE:\t AnglHyst was inicialized \n\n";
834            if((i+1)>=argc){
835                //opt1.AnglHystDefault();
836                if(opt1.debug) cout<<"Setting default value for AnglHyst option.\n\n";
837            }else{
838                
839                Int_t i0 = i+1;
840                
841                Bool_t inside = false; //Check that we go inside loop
842                
843                Int_t NPth = 0; Int_t NMA = 0; Int_t NPNO = 0;
844                opt1.AnglHystFunc(argc,argv);
845                
846            }       //if((i+1)>=argc) {...}; else
847    
848        }   //AnglHyst
849        
850        if(!strcmp(argv[i],"--DiffHyst")){
851            //opt1.opt.DH.DiffHystOn = true;
852            //opt1.DiffHystFunc(argc,argv);
853        }*/
854        if(!strcmp(argv[i],"--PamEff")){
855            //opt1.opt.PE.PamEffOn = true;
856            opt1.PamEffFunc(argc,argv);
857            /*if(opt1.debug){
858                for(Int_t i = 0; i<opt1.Hyst2D.size();i++){
859                    cout<<"\t\tName of Hystogramm - "<<opt1.Hyst2D[i].Hyst2DF.GetName()<<endl;
860                    cout<<"\t\tRoot = "<<opt1.Hyst2D[i].HI.root<<endl;
861                    cout<<"\t\tpng = "<<opt1.Hyst2D[i].HI.png<<endl;
862                    cout<<"\t\tExpType = "<<opt1.Hyst2D[i].HI.ExpType<<endl;
863                    cout<<"\t\tFile name = "<<opt1.Hyst2D[i].HI.filename<<endl;
864                }
865            }*/
866        }   //PamEff
867        
868        if(!strcmp(argv[i],"-DoTr")) opt1.opt.DoTr = true;
869        if(!strcmp(argv[i],"-v") || !strcmp(argv[i],"-verbose")){
870            if(i+1>=argc || !atoi(argv[i+1])){
871                cout<<"\n ERROR! \t -verbose needs arguments \n See --help \n";
872                if(opt1.debug) cout<<"Argument number\t"<<i<<"\n";
873                exit(1);
874            }
875            opt1.opt.verbose = atoi(argv[i+1]);
876            if(opt1.opt.verbose>3){
877                cout<<"\n ERROR! \t -verbose needs 0, 1, 2 or 3 only \n See --help \n";
878                if(opt1.debug) cout<<"Argument\t"<<opt1.opt.verbose<<"\n";
879                exit(1);
880            }
881        }
882        if(!strcmp(argv[i],"-s") || !strcmp(argv[i],"-silent")) opt1.opt.verbose = 0;
883    /*    if(!strcmp(argv[i],"-outform") || !strcmp(argv[i],"-of")){
884            Int_t i7 = i+1;
885            if(opt1.debug) cout<<"OutForm option was selected...\n\n";
886            if(i7>=argc || (strcmp(argv[i7],"jpg") && strcmp(argv[i7],"png") && strcmp(argv[i7],"bmp") && strcmp(argv[i7],"root") &&
887            strcmp(argv[i7],"txt") && strcmp(argv[i7],"all"))){
888                cout<<"\n ERROR! \t -outform needs argument \n See --help \n";
889                if(opt1.debug) cout<<"Argument number - "<<i7<<" argc = "<<argc<<"\n";
890                exit(1);
891            }
892            while(i7<argc && (!strcmp(argv[i7],"jpg") || !strcmp(argv[i7],"png") || !strcmp(argv[i7],"bmp") || !strcmp(argv[i7],"root") ||
893            !strcmp(argv[i7],"txt") || !strcmp(argv[i7],"all"))){
894                if(!strcmp(argv[i7],"jpg")) opt1.opt.ext[0] = true;
895                if(!strcmp(argv[i7],"png")) opt1.opt.ext[1] = true;
896                if(!strcmp(argv[i7],"bmp")) opt1.opt.ext[2] = true;
897                if(!strcmp(argv[i7],"root")) opt1.opt.ext[3] = true;
898                if(!strcmp(argv[i7],"txt")) opt1.opt.ext[4] = true;
899                if(!strcmp(argv[i7],"all")){
900                    opt1.opt.ext[1] = true;
901                    opt1.opt.ext[3] = true;
902                    opt1.opt.ext[4] = true;
903                }
904            i7++;
905            }
906        }*/
907        if(!strcmp(argv[i],"-ne") || (!strcmp(argv[i],"-nevents"))){
908            if(i+1>=argc || !atoi(argv[i+1])){
909                cout<<"\n ERROR! \t -nevents  needs ULong integer \n See --help \n";
910                exit(1);
911            }
912            if(atoi(argv[i+1])<0){
913                cout<<"\n ERROR! \t -nevents  needs ULong integer \n See --help \n";
914                exit(1);
915            }
916            opt1.opt.Neve = atoi(argv[i+1]);
917        }
918        if(!strcmp(argv[i],"-nv") || (!strcmp(argv[i],"-Nverbose"))){
919            if(i+1>=argc || !atoi(argv[i+1])){
920                cout<<"\n ERROR! \t -nevents  needs ULong integer \n See --help \n";
921                exit(1);
922            }
923            if(atoi(argv[i+1])<0){
924                cout<<"\n ERROR! \t -nevents  needs ULong integer \n See --help \n";
925                exit(1);
926            }
927            opt1.opt.Nverb = atoi(argv[i+1]);
928        }
929        if(!strcmp(argv[i],"-ns") || (!strcmp(argv[i],"-Nstart"))){
930            if(i+1>=argc || !atoi(argv[i+1])){
931                cout<<"\n ERROR! \t -Nstart needs ULong integer \n See --help \n";
932                exit(1);
933            }
934            if(atoi(argv[i+1])<0){
935                cout<<"\n ERROR! \t -nstart  needs ULong integer \n See --help \n";
936                exit(1);
937            }
938            opt1.opt.Nstart = atoi(argv[i+1]);
939        }
940        if(!strcmp(argv[i],"--Show")){
941            Int_t a=0;
942            if(opt1.debug) a=1;
943            if(i>5+a){
944                cout<<"\n Error!: \t --Show must be used as single option without any other \n See --help \n";
945                if(opt1.debug) cout<<"Argument number - "<<i<<" argc = "<<argc<<" arg0 - "<<argv[0]<<" arg1 - "<<argv[1]<<"\n";
946                exit(1);
947            }
948            if(opt1.debug) cout<<"\nShow Option was selected...\n\n";
949            opt1.opt.SP.Existed = true;
950            if(!strcmp(argv[i+1],"-nevents")) opt1.opt.SP.Nev = true;
951            if(!strcmp(argv[i+1],"-StartTime")) opt1.opt.SP.ST = true;
952            if(!strcmp(argv[i+1],"-EndTime")) opt1.opt.SP.ET = true;
953            if(!strcmp(argv[i+1],"-FullInformation")){
954                Int_t j6 = i+2;
955                if(opt1.debug) cout<<"\nFull Information option was selected...\n\n";
956                if(j6>=argc){
957                    cout<<"\n ERROR! \t -FullInformation needs argument \n See --help \n";
958                    exit(1);
959                }
960                while(j6<argc && (!strcmp(argv[j6],"-absTime") || !strcmp(argv[j6],"-event") || !strcmp(argv[j6],"-UTS") || !strcmp(argv[j6],"-Track"))){
961                    if(!strcmp(argv[j6],"-event")){
962                        if(j6+1>=argc || !atoi(argv[j6+1])){
963                            cout<<"\n ERROR! \t -nevents needs argument \n See --help \n";
964                            exit(1);
965                        }
966                        opt1.opt.SP.FI.NEvent = atoi(argv[j6+1]);
967                        if(opt1.debug) cout<<"Event number\t"<<opt1.opt.SP.FI.NEvent<<"\n";
968                        j6++;
969                    }
970                    if(!strcmp(argv[j6],"-absTime")){
971                        if(j6+1>=argc || !atoi(argv[j6+1])){
972                            cout<<"\n ERROR! \t -absTime needs argument \n See --help \n";
973                            exit(1);
974                        }
975                        opt1.opt.SP.FI.AbsTime = atoi(argv[j6+1]);
976                        if(opt1.debug) cout<<"Absolute time of Event\t"<<opt1.opt.SP.FI.AbsTime<<"\n";
977                        j6++;
978                    }
979                    if(!strcmp(argv[j6],"-UTS")){
980                        if(j6+1>=argc){
981                            cout<<"\n ERROR! \t -UTS needs argument \n See --help \n";
982                            exit(1);
983                        }
984                        opt1.opt.SP.FI.UTS = argv[j6+1];
985                        if(opt1.debug) cout<<"UTS of Event\t"<<opt1.opt.SP.FI.NEvent<<"\n";
986                        j6++;
987                    }
988                    if(!strcmp(argv[j6],"-Track")) opt1.opt.SP.FI.WT = true;
989                j6++;
990                }
991            }
992            if(!strcmp(argv[i+1],"-absTime")){
993                if(i+2>=argc || !atoi(argv[i+2])){
994                    cout<<"\n ERROR! \t -event needs integer argument \n See --help \n";
995                    exit(1);
996                }
997                opt1.opt.SP.AbsTime = atoi(argv[i+2]);
998                if(opt1.debug) cout<<"AbsTime for convert is "<<opt1.opt.SP.AbsTime<<"\n";
999            }
1000            if(!strcmp(argv[i+1],"-UTS")){
1001                if(i+2>=argc){
1002                    cout<<"\n ERROR! \t -UTS needs argument \n See --help \n";
1003                    exit(1);
1004                }
1005                opt1.opt.SP.UTS = argv[i+2];
1006                if(opt1.debug) cout<<"UTS for convert is "<<opt1.opt.SP.UTS<<"\n";
1007            }
1008        }
1009    }
1010    
1011    //Bool_t OptDef = false;
1012    //for(Int_t t=0; t<5; t++) if (opt1.opt.ext[t] == true) OptDef = true;
1013    //if(!OptDef) opt1.opt.ext[1] = true;
1014    /*
1015    if(opt1.debug){
1016        if(opt1.opt.SP.Existed){
1017            cout<<"\nShow options are:\n\n";
1018            cout<<"number of events - "<<opt1.opt.SP.Nev<<";\t Start time - "<<opt1.opt.SP.ST<<";\t End time - "<<opt1.opt.SP.ET<<";\n";
1019        }else{
1020            cout<<"\nOptions are:\n\n";
1021            cout<<"DoTrack - "<<opt1.opt.DoTr<<";\t verbose - "<<opt1.opt.verbose<<";\n";
1022            cout<<"output file:\t";
1023            for(Int_t t=0; t<5; t++) cout<<opt1.opt.ext[t]<<",  ";
1024            cout<<"\n\n";
1025        }
1026    }*/
1027    
1028    if(!opt1.debug) Calculate(argv[2],outdir,AnglHist,DiffHist,PamEff,DoTr,PamAngTime,PamExp,opt1);
1029    
1030  }  }
1031    
1032  //Function to Convert Rigidity to Energy  //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,
1055   return R;   return R;
1056  }  }
1057    
1058  void Calculate(char* dirname, string outdir, Bool_t AnglHist, Bool_t DiffHist, Bool_t PamEff, Bool_t DoTr, Bool_t PamAngTime, Bool_t PamExp){  
1059        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){
1060    
1061      /**********************************************/      /**********************************************/
1062      //First initialization, general for all purposes      //First initialization, general for all purposes
1063      /**********************************************/      /**********************************************/
1064            
1065      //gStyle->SetOptStat(111111);      if(opt.opt.verbose>=2) cout<<"\nGeneral variables initialisation...\n\n";
1066      TSystemDirectory *workdir = new TSystemDirectory("workdir",dirname);      
1067      TList *flist=workdir->GetListOfFiles();      TSystemDirectory* workdir = new TSystemDirectory("workdir",dirname);
1068        //TSystemDirectory* workdir2 = new TSystemDirectory("workdir2","/data01/_3/607_3/");
1069        TList* flist = workdir->GetListOfFiles();
1070      PamLevel2* pam_events = new PamLevel2();      PamLevel2* pam_events = new PamLevel2();
1071      PamelaOrientation* PO = new PamelaOrientation();      PamelaOrientation* PO = new PamelaOrientation();
1072      TTree *T = pam_events->GetPamTree(flist,"treename");      TTree *T = pam_events->GetPamTree(flist,"treename");
1073      ULong_t nevents = T->GetEntries();      ULong_t nevents = T->GetEntries();
1074    
1075        if(opt.opt.verbose>=1) cout<<"\n Nevents = "<<nevents<<"\n\n";
1076        
1077            
     cout<<"Number of events: "<<nevents<<endl;  
1078      /********************************************************************************/      /********************************************************************************/
1079      /*****NEED TO  CHANGE FOR OTHER COMPUTERS****************************************/      /*****NEED TO  CHANGE FOR OTHER COMPUTERS****************************************/
1080      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/");
1081      /********************************************************************************/      /********************************************************************************/
1082            
1083      Int_t nz = 6; Float_t zin[6];      Int_t nz = 6; Float_t zin[6];
1084      for (Int_t ip=0;ip<nz;ip++) {zin[ip] = pam_events->GetToFLevel2()->GetZTOF(pam_events->GetToFLevel2()->GetToFPlaneID(ip));cout<<zin[ip]<<endl;}      for (Int_t ip=0;ip<nz;ip++) {zin[ip] = pam_events->GetToFLevel2()->GetZTOF(pam_events->GetToFLevel2()->GetToFPlaneID(ip));cout<<zin[ip]<<endl;}
1085      Trajectory *tr = new Trajectory(nz,zin);      Trajectory *tr = new Trajectory(nz,zin);
1086        
1087        if(opt.opt.verbose>=2) cout<<"\n Other Variables initialisation...\n\n";
1088    
1089      PamTrack *track;      PamTrack *track;
1090            
1091      Int_t ntr;      Int_t ntr;                  Float_t Argv = 0;       Double_t PamAzim = 0;
1092      Float_t Argv = 0;      Double_t PamZenit = 0;      Double_t MyAzim = 0;    Double_t MyZenit = 0;
1093      Double_t PamAzim = 0;      Double_t Px = 0;            Double_t Py = 0;        Double_t Pz = 0;
1094      Double_t PamZenit = 0;      
1095      Double_t MyAzim = 0;      //Peremennye dlya opredeleniya pervogo vhoda v cycly v PamExp
1096      Double_t MyZenit = 0;      
1097      Double_t Px;      ULong_t PHE1 = 0;           ULong_t PLE1 = 0;       ULong_t PAE1 = 0;
1098      Double_t Py;      ULong_t EHE1 = 0;           ULong_t ELE1 = 0;       ULong_t EAE1 = 0;
1099      Double_t Pz;      ULong_t pHE1 = 0;           ULong_t pLE1 = 0;       ULong_t pAE1 = 0;
1100            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");  
         }  
       
     }  
1101            
1102      /*********************************************************************************/      /*********************************************************************************/
1103      //Histogramms for Count and exposition of Angles (Pitch, Pamela's Main axis, etc )      //Histogramms for Count and exposition of Angles (Pitch, Pamela's Main axis, etc )
1104      /*********************************************************************************/      /*********************************************************************************/
1105            
1106      //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);  
       
     //}  
1107            
1108      /***************************************************************************************/      /***************************************************************************************/
1109      //Histogramms for diferences in results between calculating using DoTrack2      //Histogramms for differences in results between calculating using DoTrack2
1110      //and using axv[0], ayv[0]      //and using axv[0], ayv[0]
1111      /**************************************************************************************/      /**************************************************************************************/
1112            
1113      //if(DiffHist){      if(opt.opt.verbose>=2) cout<<"\n--DiffHyst Hystogramm block...\n\n";
1114        
     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);  
1115            
1116      /***************************************************************************************/      /***************************************************************************************/
1117      //Initialization histogramms for count of protons depends from Energy. Volodia.      //Initialization histogramms for count of protons depends from Energy. Volodia.
1118      /***************************************************************************************/      /***************************************************************************************/
1119  /*          
1120      TH1F* proton10log= new TH1F("Protons Flux","",80,-1.,3.);  //    TH1F* proton10log= new TH1F("Protons Flux","",80,-1.,3.);
1121      proton10log->GetXaxis()->SetTitle("log10(Energy) , GeV");  //    proton10log->GetXaxis()->SetTitle("log10(Energy) , GeV");
1122      TH1F* proton1log= new TH1F("Protons Flux","",80,-1.,3.);  //    TH1F* proton1log= new TH1F("Protons Flux","",80,-1.,3.);
1123      TH1F* binw= new TH1F("w","",80,-1.,3.);      //TH1F* binw= new TH1F("w","",80,-1.,3.);
1124      proton1log->GetXaxis()->SetTitle("log10(Energy) , GeV");      TH1F binw;
1125      proton1log->SetLineColor(kRed);      binw.SetBins(80,-1.,3.);
1126      TCanvas *plogCanvas = new TCanvas("proton flux","protonflux", 800, 800);  //    proton1log->GetXaxis()->SetTitle("log10(Energy) , GeV");
1127    //    proton1log->SetLineColor(kRed);
1128    //    TCanvas *plogCanvas = new TCanvas("proton flux","protonflux", 800, 800);
1129      // this is bins wide to calculate Flux      // this is bins wide to calculate Flux
1130      Float_t x;      //Float_t x;
1131      Float_t xmin, xmax;      //Float_t xmin, xmax;
1132      Float_t binwide[80];      //Float_t binwide[80];
1133      for (Int_t l=0 ; l<80; l++) {      for (Int_t l=0 ; l<80; l++) {
1134          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.));
1135            //binwide[l]=-1.+0.05*l,pow(10.,0.05*(Float_t)(l+1)-1.)- pow(10.,0.05*(Float_t)(l)-1.)
1136      }      }
1137        
1138      for (Int_t l=0 ; l<80; l++) {      //for (Int_t l=0 ; l<80; l++) {
1139      binwide[l]=binw->GetBinContent(l+1);      //binwide[l]=binw->GetBinContent(l+1);
1140      }  //    }
1141  */      
1142      /****************************************************************************/      /****************************************************************************/
1143      //Geometrical Factor. Volodia      //Geometrical Factor. Volodia
1144      /****************************************************************************/      /****************************************************************************/
# Line 305  void Calculate(char* dirname, string out Line 1161  void Calculate(char* dirname, string out
1161      //General loop      //General loop
1162      /****************************************************************************/      /****************************************************************************/
1163    
1164      for(ULong_t i=0; i<nevents; i++){      if(opt.opt.Neve == 0) opt.opt.Neve = nevents;
1165    
1166        cout<<opt.opt.Neve<<endl;
1167    
1168        for(ULong_t i=opt.opt.Nstart; i<opt.opt.Neve; i++){
1169            
1170            //cout<<i<<endl;
1171          T->GetEntry(i);          T->GetEntry(i);
1172            Bool_t OneTrack = false;
1173    
1174            if(pam_events->GetTrkLevel2()->GetNTracks()==1) OneTrack = true;
1175                    
1176          if(pam_events->GetTrkLevel2()->GetNTracks()==1){          Int_t M0DE = pam_events->GetOrbitalInfo()->mode;  //0 is zero here
1177                        Float_t Bx = -pam_events->GetOrbitalInfo()->Bdown;                      //don't need for PamExp ExpOnly for all geography areas
1178              /*************************************************************************/          Float_t By = pam_events->GetOrbitalInfo()->Beast;                       //don't need for PamExp ExpOnly for all geography areas
1179              //Getting general parameters          Float_t Bz = pam_events->GetOrbitalInfo()->Bnorth;                      //don't need for PamExp ExpOnly for all geography areas
1180              /*************************************************************************/          Float_t Babs = pam_events->GetOrbitalInfo()->Babs;                      //don't need for PamExp ExpOnly for all geography areas
1181                        Float_t L = pam_events->GetOrbitalInfo()->L;                            //don't need for PamExp ExpOnly for all geography areas
1182              track = pam_events->GetTrack(0);          Float_t Alt = pam_events->GetOrbitalInfo()->alt;
1183              Int_t M0DE = pam_events->GetOrbitalInfo()->mode;  //0 is zero here          Float_t Lat = pam_events->GetOrbitalInfo()->lat;
1184              Float_t Bx = -pam_events->GetOrbitalInfo()->Bdown;          Float_t Lon = pam_events->GetOrbitalInfo()->lon;
1185              Float_t By = pam_events->GetOrbitalInfo()->Beast;          Float_t LieveTime = 0.16*pam_events->GetTrigLevel2()->dltime[0];
1186              Float_t Bz = pam_events->GetOrbitalInfo()->Bnorth;          Float_t DeathTime = 0.01*pam_events->GetTrigLevel2()->dltime[1];
1187              Float_t Babs = pam_events->GetOrbitalInfo()->Babs;          Float_t rigev = 0;
1188              Float_t L = pam_events->GetOrbitalInfo()->L;          Float_t detr = 0;
1189              Float_t rigev = 1/track->GetTrkTrack()->GetDeflection();          Int_t hm0; Int_t hm1;
1190                        Int_t hs0; Int_t hs1;
1191              /*************************************************************************/  
1192              //Track selection          ntr = 0;    
1193              /*************************************************************************/  
1194                        Double_t A1; Double_t A2; Double_t A3;
1195              ntr = 0;          Double_t SB = 1000; Double_t SA; Double_t SC; Double_t ZC;
1196              if ((track->GetTrkTrack()->chi2<=0) || (track->GetTrkTrack()->chi2>=500)) ntr=-1;  
1197              for(Int_t ii=1;ii<=12;ii++){          if(OneTrack){
1198                  if ((track->GetToFTrack()->beta[ii]<0) || (track->GetToFTrack()->beta[ii]>99)) ntr=-1;              track = pam_events->GetTrack(0);                            //don't need for PamExp option
1199              }              rigev = 1/track->GetTrkTrack()->GetDeflection();            //dont need for AnglHyst, DiffHyst and PamExp options;
1200              if((track->GetTrkTrack()->GetNX()<=4)&&(track->GetTrkTrack()->GetNY()<4)) ntr=-1;              detr = track->GetTrkTrack()->GetDEDX();
1201              if((M0DE!=0)&&(M0DE!=1)&&(M0DE!=6)&&(M0DE!=2)&&(M0DE!=3)&&(M0DE!=8)&&(M0DE!=4)) ntr=-1;              hm0 = pam_events->GetAcLevel2()->hitmap[0];
1202              if(sqrt(pow(Bx,2)+pow(By,2)+pow(Bz,2))<0.001) ntr=-1;              hm1 = pam_events->GetAcLevel2()->hitmap[1];
1203              //if(rigev<0.600) ntr=-1;              hs0 = pam_events->GetAcLevel2()->hitstatus[0];
1204                            hs1 = pam_events->GetAcLevel2()->hitstatus[1];
1205              if(ntr!=-1 || PamExp){              
1206                            if(!OneTrack) ntr = -1; else{
1207                  /*******************************************************************************/                  //cout<<"OneTrack"<<endl;
1208                  //Transit reference system                  Float_t ichi2lim = fabs(track->GetTrkTrack()->GetDeflection())*1.85+3.0;
1209                  /*******************************************************************************/                  if ((track->GetTrkTrack()->chi2<=0) || (track->GetTrkTrack()->chi2>=pow(ichi2lim,4.))) ntr=-1;
1210                                /*for(Int_t ii=1;ii<=12;ii++){
1211                  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;
1212                  TMatrixD Dij = PO->GreenwichtoGEO(pam_events->GetOrbitalInfo()->lat,pam_events->GetOrbitalInfo()->lon,Fij);                  }*/
1213                  TMatrixD Iij = PO->ColPermutation(Dij);                  if ((track->GetToFTrack()->beta[12]<0) || (track->GetToFTrack()->beta[12]>99)) ntr=-1;
1214                                if((track->GetTrkTrack()->GetNX()<=4)&&(track->GetTrkTrack()->GetNY()<4)) ntr=-1;
1215                  /*****************************************************************************/                  if(!track->GetTrkTrack()->IsInsideCavity()) ntr=-1;
1216                  //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;
1217                  //using axv[0] and ayv[0] variables                  if((M0DE!=0)&&(M0DE!=1)&&(M0DE!=6)&&(M0DE!=2)&&(M0DE!=3)&&(M0DE!=8)&&(M0DE!=4)) ntr=-1;
1218                  /*****************************************************************************/                  }
1219                            
1220                  if(!DoTr || DiffHist){              if(ntr!=-1){
1221                            
1222                if(!opt.opt.DoTr){
1223                    
1224                  Double_t Aaxv = TMath::Abs(track->GetTrkTrack()->axv[0])*TMath::DegToRad();                  Double_t Aaxv = TMath::Abs(track->GetTrkTrack()->axv[0])*TMath::DegToRad();
1225                  Double_t Aayv = TMath::Abs(track->GetTrkTrack()->ayv[0])*TMath::DegToRad();                  Double_t Aayv = TMath::Abs(track->GetTrkTrack()->ayv[0])*TMath::DegToRad();
1226                  PamZenit = TMath::RadToDeg()*asin(sqrt(pow(sin(Aayv), 2) + pow(sin(Aaxv), 2)));                  PamZenit = TMath::RadToDeg()*asin(sqrt(pow(sin(Aayv), 2) + pow(sin(Aaxv), 2)));
1227                    
1228                  Double_t axv = -track->GetTrkTrack()->axv[0] * TMath::DegToRad();                  Double_t axv = -track->GetTrkTrack()->axv[0] * TMath::DegToRad();
1229                  Double_t ayv = -track->GetTrkTrack()->ayv[0] * TMath::DegToRad();                  Double_t ayv = -track->GetTrkTrack()->ayv[0] * TMath::DegToRad();
1230                  Double_t angle = atan(sin(TMath::Abs(ayv))/sin(TMath::Abs(axv))) * TMath::RadToDeg();                  Double_t angle = atan(sin(TMath::Abs(ayv))/sin(TMath::Abs(axv))) * TMath::RadToDeg();
1231                    
1232                  PamAzim =  360. - angle;                  PamAzim =  360. - angle;
1233                  if(axv>=0 && ayv >=0) PamAzim = angle;                  if(axv>=0 && ayv >=0) PamAzim = angle;
1234                  if(axv<0 && ayv >0) PamAzim = 180. - angle;                  if(axv<0 && ayv >0) PamAzim = 180. - angle;
1235                  if(axv<0 && ayv <0) PamAzim = 180. + angle;                  if(axv<0 && ayv <0) PamAzim = 180. + angle;
1236                    
1237                  PamAzim = PamAzim * TMath::DegToRad();                  PamAzim = PamAzim * TMath::DegToRad();
1238                  PamZenit = (180 - PamZenit) * TMath::DegToRad();                  PamZenit = (180 - PamZenit) * TMath::DegToRad();
1239                  Px = sin(axv);                  Px = sin(axv);
1240                  Py = sin(ayv);                  Py = sin(ayv);
1241                  Pz = cos(PamZenit);                  Pz = cos(PamZenit);
1242                                    
1243                    //cout<<"PamAzimuth = "<<PamAzim<<" PamZenit = "<<PamZenit<<endl;
1244                                    
1245                  }              }
1246                            
1247                  /*****************************************************************************/              /*****************************************************************************/
1248                  //Calculate Zenit and Azimutal angle and vector of particle flight in Pamela              //Calculate Zenit and Azimutal angle and vector of particle flight in Pamela
1249                  //using DoTrack2 procedure              //using DoTrack2 procedure
1250                  /*****************************************************************************/              /*****************************************************************************/
1251                            
1252                  if(DoTr){              if(opt.opt.DoTr){
1253                            
1254                  track->GetTrkTrack()->DoTrack2(tr);                  track->GetTrkTrack()->DoTrack2(tr);
1255                  Double_t E11x = tr->x[0];                  Double_t E11x = tr->x[0];
1256                  Double_t E11y = tr->y[0];                  Double_t E11y = tr->y[0];
1257                  Double_t E11z = zin[0];                  Double_t E11z = zin[0];
1258                  Double_t E22x = tr->x[3];                  Double_t E22x = tr->x[3];
1259                  Double_t E22y = tr->y[3];                  Double_t E22y = tr->y[3];
1260                  Double_t E22z = zin[3];                  Double_t E22z = zin[3];
1261                  Double_t norm = sqrt(pow(E22x-E11x,2)+pow(E22y-E11y,2)+pow(E22z-E11z,2));                  Double_t norm = sqrt(pow(E22x-E11x,2)+pow(E22y-E11y,2)+pow(E22z-E11z,2));
1262                  MyAzim = TMath::RadToDeg()*atan(TMath::Abs(E22y-E11y)/TMath::Abs(E22x-E11x));                  MyAzim = TMath::RadToDeg()*atan(TMath::Abs(E22y-E11y)/TMath::Abs(E22x-E11x));
1263                  MyAzim = MyAzim;//-180;                  MyAzim = MyAzim;//-180;
1264                  MyZenit = TMath::RadToDeg()*acos((E22z-E11z)/norm);                  MyZenit = TMath::RadToDeg()*acos((E22z-E11z)/norm);
1265                  if(E22x-E11x>=0 && E22y-E11y <0) MyAzim =  360. - MyAzim;                  if(E22x-E11x>=0 && E22y-E11y <0) MyAzim =  360. - MyAzim;
1266                  if(E22x-E11x>=0 && E22y-E11y >=0) MyAzim = MyAzim;                  if(E22x-E11x>=0 && E22y-E11y >=0) MyAzim = MyAzim;
1267                  if(E22x-E11x<0 && E22y-E11y >0) MyAzim = 180. - MyAzim;                  if(E22x-E11x<0 && E22y-E11y >0) MyAzim = 180. - MyAzim;
1268                  if(E22x-E11x<0 && E22y-E11y <0) MyAzim = 180. + MyAzim;                  if(E22x-E11x<0 && E22y-E11y <0) MyAzim = 180. + MyAzim;
1269                  Px = (E22x-E11x)/norm;                  Px = (E22x-E11x)/norm;
1270                  Py = (E22y-E11y)/norm;                  Py = (E22y-E11y)/norm;
1271                  Pz = (E22z-E11z)/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  
1272                                    
1273                  /******************************************************************************/              }
1274                  //Proton selection. Volodia.                  
1275                  /******************************************************************************/              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);
1276  /*                            TMatrixD Dij = PO->GreenwichtoGEO(pam_events->GetOrbitalInfo()->lat,pam_events->GetOrbitalInfo()->lon,Fij);
1277                  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){  
1278                                    
1279                  if(L>6){              TMatrixD Eij = PO->PamelatoGEO(Iij,Px,Py,Pz);               //Don't necessary for PamExp
1280                  Float_t Mass = 0;                  
1281                  Int_t detr = 0;              A1 = Iij(0,2);
1282                  Float_t Ekin = 0;              A2 = Iij(1,2);
1283                  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);
1284                      if(14.9/L/L > 0. && 14.9/L/L < 1. & Babs >0.24) {                    
1285                          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
1286                              Ekin=ConvertR2T(rigev,0.511,-1.);              SA = PO->GetPitchAngle(A1,A2,A3,Bx,By,Bz);                  // Angle between Pamela's main axiz and B
1287                              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
1288                              if (Ekin<1) PamAngEfflow->Fill(TMath::RadToDeg()*PamAzim,TMath::RadToDeg()*PamZenit);              
1289                          }//electron2->Fill(log10(fabs(rigev)));              }
1290                                    }
1291                          if(detr  < (15*exp(-1.6*(fabs(rigev)-0.5))+4)  && rigev >0.5  || (detr < 30 && fabs(rigev) <0.5)) {          
1292                              Ekin=ConvertR2T(rigev,0.938,-1.);          //cout<<"OneTrack"<<endl;
1293                              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{
1294                              if (Ekin<1) PamAngEfflow->Fill(TMath::RadToDeg()*PamAzim,TMath::RadToDeg()*PamZenit);              //cout<<"US.size = "<<opt.US.size()<<endl;
1295                for(Int_t jz = 0; jz<opt.PS.size();jz++){
1296                    for(Int_t jy = 0; jy<opt.US.size();jy++){
1297                        for(Int_t jx = 0; jx<opt.US[jy].USP.size();jx++){
1298                            //if(SB<400) cout<<"Pitch = "<<SB<<" Babs = "<<Babs<<" L = "<<L<<" Alt = "<<Alt<<endl;
1299                            //cout<<"Fill Hyst1D["<<thg<<"]"<<endl;
1300                            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){
1301                                for(Int_t hg = 0; hg<opt.US[jy].NAH1D.size(); hg++){
1302                                    Int_t thg = opt.US[jy].NAH1D[hg];
1303                                    for(Int_t ps = 0; ps<opt.PS[jz].NAH1D.size();ps++){
1304                                        if(thg==opt.PS[jz].NAH1D[ps]){
1305                                            opt.Hyst1D[thg].HI.LieveTime+= 0.16*pam_events->GetTrigLevel2()->dltime[0];
1306                                            opt.Hyst1D[thg].HI.DeathTime+= 0.01*pam_events->GetTrigLevel2()->dltime[1];
1307                                            //cout<<"Pitch = "<<SB<<" Babs = "<<Babs<<" L = "<<L<<" Alt = "<<Alt<<endl;
1308                                            //cout<<"Fill Hyst1D["<<thg<<"]"<<endl;
1309                                        }
1310                                    }
1311                                }
1312                                for(Int_t hg = 0; hg<opt.US[jy].NAH2D.size(); hg++){
1313                                    Int_t thg = opt.US[jy].NAH2D[hg];
1314                                    for(Int_t ps = 0; ps<opt.PS[jz].NAH2D.size();ps++){
1315                                        if(thg==opt.PS[jz].NAH2D[ps]){
1316                                            opt.Hyst2D[thg].HI.LieveTime+= 0.16*pam_events->GetTrigLevel2()->dltime[0];
1317                                            opt.Hyst2D[thg].HI.DeathTime+= 0.01*pam_events->GetTrigLevel2()->dltime[1];
1318                                        }
1319                                    }
1320                                }
1321                                for(Int_t hg = 0; hg<opt.US[jy].NAH3D.size(); hg++){
1322                                    Int_t thg = opt.US[jy].NAH3D[hg];
1323                                    for(Int_t ps = 0; ps<opt.PS[jz].NAH3D.size();ps++){
1324                                        if(thg==opt.PS[jz].NAH3D[ps]){
1325                                            opt.Hyst3D[thg].HI.LieveTime+= 0.16*pam_events->GetTrigLevel2()->dltime[0];
1326                                            opt.Hyst3D[thg].HI.DeathTime+= 0.01*pam_events->GetTrigLevel2()->dltime[1];
1327                                        }
1328                                    }
1329                                }
1330                                jx = opt.US[jy].USP.size();
1331                          }                          }
1332                      }                      }
1333                  }                  }
1334                                }
1335                  }          }
1336                            //cout<<ntr<<endl;
1337                  /***********************************************************************/          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){  
1338    
1339                  /***********************************************************************/              if(ntr!=-1){
1340                  //Filling angular histogramms for Brazil              
1341                  /***********************************************************************/                  Bool_t verb = false;
1342                    if(opt.opt.Nverb!=-1) if(i>=opt.opt.Nverb) verb = true;
1343                  if((Babs<0.19)&&(L<1.75)){                  //if(opt.opt.Nverb==-1) verb = false;
1344                          Argv = 0.16*pam_events->GetTrigLevel2()->dltime[0];  
1345                          PitchExpositionBrazil->Fill(SB,Argv);                  if(verb) cout<<"\n\ni = "<<i<<"\n\n"<<endl;
1346                          MainAxisPamelaExpositionBrazil->Fill(SA,Argv);                                          
1347                          PitchExpositionNoOrientationBrazil->Fill(SC,Argv);                  /*****************************************************************************/
1348                          PitchCountBrazil->Fill(SB);                  // CALO Selection
1349                          MainAxisPamelaCountBrazil->Fill(SA);                  /*****************************************************************************/
1350                          PitchCountNoOrientationBrazil->Fill(SC);                  
1351                          ZenitAngle->Fill(ZC);                  Float_t caloqpre=track->GetCaloTrack()->qpre;
1352                          PitchvsLatInBrazil->Fill(pam_events->GetOrbitalInfo()->lat,SA);                  Int_t   calonpre=track->GetCaloTrack()->npre;
1353                          LieveTimevsLatInBrazil->Fill(pam_events->GetOrbitalInfo()->lat,Argv);                  Float_t caloqtrack=track->GetCaloTrack()->qtrack;
1354                  }                  Float_t caloqmax=pam_events->GetCaloLevel2()->qmax;
1355                                    Int_t   dncore=track->GetCaloTrack()->ncore;
1356                  /***********************************************************************/                  Int_t   lepton = 0;
1357                  //Filling angulars histogramms for areas without Brazil                  Bool_t  ACSelection = false;
1358                  /***********************************************************************/                  
1359                    if((hm0 & hs0)==0 && (hm1 & hs1)==0) ACSelection = true;
1360                    //cout<<"hm0 = "<<hm0<<" hm1 = "<<hm1<<" hs0 = "<<" hs1 = "<<hs1<<endl;
1361                    
1362                    Float_t betaev = 0.25*(track->GetToFTrack()->beta[0]+track->GetToFTrack()->beta[1]+track->GetToFTrack()->beta[2]+track->GetToFTrack()->beta[3]);
1363                    Float_t Ekin = 0;
1364    
1365                    if(detr<2. && Babs>0.24 && betaev>0.2){
1366                        if(dncore > (pow(fabs(rigev),0.4)*700.-400.) && fabs(rigev)<11.2) lepton = 1;
1367                        if(dncore > (pow(fabs(rigev),0.38)*700.-200.) && fabs(rigev)<48.) lepton = 1;
1368                        if(dncore > 3200. && fabs(rigev)>=48.) lepton = 1;
1369                        if(caloqtrack < (115.*fabs(rigev)-50.)) lepton = 0;
1370                        if(dncore < (700.*pow(fabs(rigev),0.33)-400.) && fabs(rigev)<=70. && fabs(rigev)>32.) lepton = 0;
1371                        if(fabs(rigev)>70. && dncore < 3200.) lepton = 0;
1372                        if(dncore < (pow(fabs(rigev),0.38)*700.-200.) && fabs(rigev)>18. && fabs(rigev)<32.) lepton = 0;
1373                        if(dncore < (700.*pow(fabs(rigev),0.4)-400.) && fabs(rigev)<18.) lepton = 0;
1374                        if(fabs(rigev)<=40. && calonpre < pow(fabs(rigev),0.6)) lepton = 0;
1375                        if(fabs(rigev)<=30. && calonpre > (70.-pow(fabs(rigev),0.4))) lepton = 0;
1376                        if(fabs(rigev)>40. && calonpre < 10.) lepton = 0;
1377                        if(fabs(rigev)>30. && calonpre > 32.) lepton = 0;
1378                        if(caloqpre/caloqtrack > (exp(-fabs(rigev)/5.)+exp(-fabs(rigev)/0.1)+0.05)) lepton = 0;
1379                        if(caloqmax/caloqtrack > (0.035+0.1/fabs(rigev))) lepton = 0;
1380                    }
1381    
1382                    string Gname = "";
1383                    //if(verb)for(Int_t o=0;o<opt.Hyst1D.size(); o++) cout<<"Name of Hyst1D["<<o<<"] is "<<opt.Hyst1D[o].HI.filename<<endl;
1384                    //if(verb)for(Int_t o=0;o<opt.Hyst2D.size(); o++) cout<<"Name of Hyst2D["<<o<<"] is "<<opt.Hyst2D[o].HI.filename<<endl;
1385                                    
1386                  if(Babs>0.24){                  if(verb){
1387                          Argv = 0.16*pam_events->GetTrigLevel2()->dltime[0];                  //if(rigev>1.7 && rigev<= 1.8 && L>6 && Lat>60 && Lat<-60){
1388                          PitchExpositionWithoutBrazil->Fill(SB,Argv);                      cout<<"rigev = "<<rigev<<"\tdetr = "<<detr<<endl;
1389                          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);  
                         }  
1390                  }                  }
1391                    /*cout<<"Write Rigidity, please:"<<endl;
1392                    cin>>rigev;
1393                    cout<<"Write detr, please:"<<endl;
1394                    cin>>detr;
1395                    cout<<"Write lepton, please:"<<endl;
1396                    cin>>lepton;
1397                    cout<<"Write Anticounter, please:"<<endl;
1398                    cin>>ACSelection;
1399                    cout<<"Write Babs, please:"<<endl;
1400                    cin>>Babs;
1401                    cout<<"Write L, please:"<<endl;
1402                    cin>>L;
1403                    cout<<"Write Altitude, please:"<<endl;
1404                    cin>>Alt;*/
1405                                    
1406                  } // if(AnglHist)                  Double_t Rmin = 0;
1407                                Double_t Rmax = 0;
1408              } // if(PamEff||AnglHist)                  Double_t Detrmin = 0;
1409                                Double_t Detrmax = 0;
1410                    for(Int_t k=0;k<opt.SS.size();k++){
1411                        //Int_t NSS = opt.Hyst1D[k].NSS;
1412                        Int_t ntr = 0;
1413                        Int_t u = 0;
1414                        while(u<opt.SS[k].SP.size()){
1415                            Double_t Detrmin = 0;
1416                            Double_t Detrmax = 0;
1417                            //if(verb)cout<<"\t\trigev = "<<rigev<<endl;
1418                            //if(verb)cout<<"\t\tSelectionRmin["<<u<<"] = "<<opt.SS[k].SP[u].Rmin<<"\tSelectionRmax["<<u<<"] = "<<opt.SS[k].SP[u].Rmax<<endl;
1419                            if(rigev>opt.SS[k].SP[u].Rmin && rigev<=opt.SS[k].SP[u].Rmax){
1420                                //if(verb)cout<<"\t\tdetrminf["<<u<<"] = "<<opt.SS[k].SP[u].detrminf<<"\tdetrmaxf["<<u<<"] = "<<opt.SS[k].SP[u].detrmaxf<<endl;
1421                                if(opt.SS[k].SP[u].detrminf!=0){
1422                                    switch ( opt.SS[k].SP[u].detrminf ) {
1423                                        case 1:{Detrmin = 1/pow(fabs(rigev)+0.46,5)+1.1;if(verb){cout<<"\t\tCase 1: Detrmin = "<<Detrmin<<endl;}break;}
1424                                        case 2:{Detrmin = 9.0*exp(-5.6*(fabs(rigev)-0.385))+3.7;if(verb){cout<<"\t\tCase 2: Detrmin = "<<Detrmin<<endl;}break;}
1425                                        case 3:{Detrmin = 15*exp(-1.6*(fabs(rigev)-0.5))+4.;if(verb){cout<<"\t\tCase 3: Detrmin = "<<Detrmin<<endl;}break;}
1426                                    };
1427                                }else {Detrmin = opt.SS[k].SP[u].detrmore;if(verb){cout<<"\t\tDetrmore: Detrmin = "<<Detrmin<<endl;}}
1428                                if(opt.SS[k].SP[u].detrmaxf!=0){
1429                                    switch ( opt.SS[k].SP[u].detrmaxf ) {
1430                                        case 1:{Detrmax = 1/pow(fabs(rigev)+0.46,5)+1.1;if(verb){cout<<"\t\tCase 1: Detrmax = "<<Detrmax<<endl;}break;}
1431                                        case 2:{Detrmax = 9.0*exp(-5.6*(fabs(rigev)-0.385))+3.7;if(verb){cout<<"\t\tCase 2: Detrmax = "<<Detrmax<<endl;}break;}
1432                                        case 3:{Detrmax = 15*exp(-1.6*(fabs(rigev)-0.5))+4.;if(verb){cout<<"\t\tCase 3: Detrmax = "<<Detrmax<<endl;}break;}
1433                                    };
1434                                }else {Detrmax = opt.SS[k].SP[u].detrless;if(verb){cout<<"\t\tDetrless: Detrmax = "<<Detrmax<<endl;}}
1435                                Int_t Lept;
1436                                //if(verb)cout<<"\t\t Part Type is "<<opt.SS[k].PartType<<endl;
1437                                if(opt.SS[k].PartType==1 || opt.SS[k].PartType==4) Lept =1; else Lept = 0;
1438                                if(verb)cout<<"\t\tlepton = "<<lepton<<"\tLept = "<<Lept<<endl;
1439                                if(!(Detrmax==0 && Detrmin==0) && (detr<=Detrmin || detr>Detrmax)) ntr = -1;
1440                                //if(verb)cout<<"\t\tdetr = "<<detr<<"\tDetrmin = "<<Detrmin<<"\tDetrmax = "<<Detrmax<<endl;
1441                                if(opt.SS[k].SP[u].Calo && lepton==Lept) ntr=-1;
1442                                if(opt.SS[k].SP[u].AC && !ACSelection) ntr=-1;
1443                                //if(verb)cout<<"\t\tntr = "<<ntr<<endl;
1444                            }
1445                            u++;
1446                        }//while SP  
1447                        if(ntr!=-1) {
1448                            vector<Int_t> N1D(0);
1449                            vector<Int_t> N2D(0);
1450                            Double_t Bmore = 0;
1451                            Double_t Bless = 0;
1452                            Double_t Lmore = 0;
1453                            Double_t Lless = 0;
1454                            Double_t altmore = 0;
1455                            Double_t altless = 0;
1456                            Double_t Rmore = 0;
1457                            Double_t Rless = 0;
1458                            Double_t Pmore = 0;
1459                            Double_t Pless = 0;
1460                            for(Int_t cu=0;cu<opt.SS[k].NAH1D.size();cu++){
1461                                Int_t nus=opt.Hyst1D[opt.SS[k].NAH1D[cu]].HI.NUS;
1462                                Int_t nps=opt.Hyst1D[opt.SS[k].NAH1D[cu]].HI.NPT;
1463                                for(Int_t jx = 0; jx<opt.US[nus].USP.size(); jx++){
1464                                    Bmore = opt.US[nus].USP[jx].Bmore;
1465                                    Bless = opt.US[nus].USP[jx].Bless;
1466                                    Lmore = opt.US[nus].USP[jx].Lmore;
1467                                    Lless = opt.US[nus].USP[jx].Lless;
1468                                    altmore = opt.US[nus].USP[jx].altmore;
1469                                    altless = opt.US[nus].USP[jx].altless;
1470                                    Rmore = opt.US[nus].USP[jx].Rmore;
1471                                    Rless = opt.US[nus].USP[jx].Rless;
1472                                    Pmore = opt.PS[nps].PitchMin;
1473                                    Pless = opt.PS[nps].PitchMax;
1474                                    //cout<<"Pmore = "<<Pmore<<" Pless = "<<Pless<<endl;
1475                                    if(Babs>Bmore && Babs<=Bless && L>Lmore && L<=Lless && Alt>altmore && Alt<=altless && rigev>Rmore && rigev<=Rless && SB>Pmore && SB<=Pless){
1476                                        //cout<<log10(rigev)<<endl;
1477                                        N1D.resize(N1D.size()+1);
1478                                        N1D[N1D.size()-1] = opt.SS[k].NAH1D[cu];
1479                                        if(verb){
1480                                            cout<<"nus = "<<nus<<endl;
1481                                            cout<<"Bmore = "<<Bmore<<"\tBless = "<<Bless<<endl;
1482                                            cout<<"Lmore = "<<Lmore<<"\tLless = "<<Lless<<endl;
1483                                            cout<<"N1D["<<N1D.size()-1<<"] = "<<opt.SS[k].NAH1D[cu]<<endl;
1484                                            cout<<"HystName "<<opt.Hyst1D[N1D[N1D.size()-1]].HI.filename<<endl;
1485                                            //Int_t we;
1486                                            //cin>>we;
1487                                        }
1488                                    jx = opt.US[nus].USP.size();
1489                                    }
1490                                }
1491                            }
1492                            for(Int_t cu=0;cu<opt.SS[k].NAH2D.size();cu++){
1493                                Int_t nus=opt.Hyst2D[opt.SS[k].NAH2D[cu]].HI.NUS;
1494                                Int_t nps=opt.Hyst2D[opt.SS[k].NAH2D[cu]].HI.NPT;
1495                                for(Int_t jx = 0; jx<opt.US[nus].USP.size(); jx++){
1496                                    Bmore = opt.US[nus].USP[jx].Bmore;
1497                                    Bless = opt.US[nus].USP[jx].Bless;
1498                                    Lmore = opt.US[nus].USP[jx].Lmore;
1499                                    Lless = opt.US[nus].USP[jx].Lless;
1500                                    altmore = opt.US[nus].USP[jx].altmore;
1501                                    altless = opt.US[nus].USP[jx].altless;
1502                                    Rmore = opt.US[nus].USP[jx].Rmore;
1503                                    Rless = opt.US[nus].USP[jx].Rless;
1504                                    Pmore = opt.PS[nps].PitchMin;
1505                                    Pless = opt.PS[nps].PitchMax;
1506                                    if(Babs>Bmore && Babs<Bless && L>Lmore && L<Lless && Alt>altmore && Alt<altless && rigev>Rmore && rigev<Rless && SB>Pmore && SB<Pless){
1507                                        N2D.resize(N2D.size()+1);
1508                                        N2D[N2D.size()-1] = opt.SS[k].NAH2D[cu];
1509                                        if(verb){
1510                                        //if(rigev>1.7 && rigev<= 1.8){
1511                                            cout<<"nus = "<<nus<<endl;
1512                                            cout<<"Babs = "<<Babs<<" Bmore = "<<Bmore<<"\tBless = "<<Bless<<endl;
1513                                            cout<<"L = "<<L<<" Lmore = "<<Lmore<<"\tLless = "<<Lless<<endl;
1514                                            cout<<"R = "<<rigev<<" Rmore = "<<Rmore<<"\tRless = "<<Rless<<endl;
1515                                            cout<<"SB = "<<TMath::RadToDeg()*SB<<" Pmore = "<<Pmore<<"\tPless = "<<Pless<<endl;
1516                                            cout<<"N2D["<<N2D.size()-1<<"] = "<<opt.SS[k].NAH2D[cu]<<endl;
1517                                            cout<<"HystName "<<opt.Hyst2D[N2D[N2D.size()-1]].HI.filename<<endl;
1518                                            //Int_t we;
1519                                            //cin>>we;
1520                                        }
1521                                    jx = opt.US[nus].USP.size();
1522                                    }
1523                                }
1524                            }
1525                            for(Int_t d=0;d<N1D.size();d++){
1526                                Rmin = opt.ED[opt.Hyst1D[N1D[d]].HI.NED].Emin;
1527                                Rmax = opt.ED[opt.Hyst1D[N1D[d]].HI.NED].Emax;
1528                                Detrmin = opt.ED[opt.Hyst1D[N1D[d]].HI.NED].detrmin;
1529                                Detrmax = opt.ED[opt.Hyst1D[N1D[d]].HI.NED].detrmax;
1530                                Double_t Argv = 0;
1531                                if(rigev>Rmin && rigev<=Rmax && detr>Detrmin && detr<=Detrmax){
1532                                    Int_t ExpType = opt.Hyst1D[N1D[d]].HI.ExpType;
1533                                    Double_t LTI = 0;
1534                                    Double_t DTI = 0;
1535                                    if(ExpType == 1 || ExpType == 3) LTI = 0.16*pam_events->GetTrigLevel2()->dltime[0];
1536                                    if(ExpType == 2 || ExpType == 3) DTI = 0.01*pam_events->GetTrigLevel2()->dltime[1];
1537                                    switch  ( ExpType ){
1538                                        case 1:{Argv = LTI;break;}
1539                                        case 2:{Argv = DTI;break;}
1540                                        case 3:{Argv = LTI+DTI;break;}
1541                                        case 4:{
1542                                                Argv = opt.Hyst1D[N1D[d]].HI.LieveTime;
1543                                                opt.Hyst1D[N1D[d]].HI.LieveTime = 0;
1544                                                break;
1545                                                }
1546                                        case 5:{
1547                                                Argv = opt.Hyst1D[N1D[d]].HI.DeathTime;
1548                                                opt.Hyst1D[N1D[d]].HI.DeathTime = 0;
1549                                                break;
1550                                                }
1551                                        case 6:{
1552                                                Argv = opt.Hyst1D[N1D[d]].HI.LieveTime+opt.Hyst1D[N1D[d]].HI.DeathTime;
1553                                                opt.Hyst1D[N1D[d]].HI.LieveTime = 0;
1554                                                opt.Hyst1D[N1D[d]].HI.DeathTime = 0;
1555                                                break;
1556                                                }
1557                                    };
1558                                    Double_t Xfill = 0;
1559                                    //cout<<"XType = "<<opt.Hyst1D[d].XType<<endl;
1560                                    //cout<<"Babs = "<<Babs<<endl;
1561                                    switch ( opt.Hyst1D[N1D[d]].XType ){
1562                                        case 1:{ Xfill = TMath::RadToDeg()*PamAzim;break; }
1563                                        case 2:{ Xfill = TMath::RadToDeg()*PamZenit;break; }
1564                                        case 3:{ Xfill = SB; break; }
1565                                        case 4:{ Xfill = detr; break; }
1566                                        case 5:{ Xfill = rigev; break; }
1567                                        case 6:{ Xfill = Lat; break; }
1568                                        case 7:{ Xfill = Lon; break; }
1569                                        case 8:{ Xfill = L; break; }
1570                                        case 9:{ Xfill = Babs; break; }
1571                                        case 10:{Xfill = Alt; break; }
1572                                        case 11:{Xfill = LieveTime; break; }
1573                                        case 12:{Xfill = DeathTime; break; }
1574                                        case 13:{Xfill = SC; break; }
1575                                        case 14:{Xfill = SA; break; }
1576                                        case 17:{Xfill = log10(rigev); break;}
1577                                    };
1578                                    if(verb){
1579                                        cout<<"rigev = "<<rigev<<"\tdetr = "<<detr<<endl;
1580                                        cout<<"Babs = "<<Babs<<"\tL = "<<L<<endl;
1581                                        cout<<"I Fill "<<opt.Hyst1D[N1D[d]].HI.filename<<" hystogramm"<<endl;
1582                                        cout<<"Xfill = "<<Xfill<<" d = "<<d<<endl;
1583                                    }
1584                                    opt.Hyst1D[N1D[d]].Hyst1DF.Fill(Xfill);
1585                                    if(ExpType!=0) opt.Hyst1D[N1D[d]].Hyst1DFExp.Fill(Xfill,Argv);
1586                                    if(verb){
1587                                        Int_t intt;
1588                                        cin>>intt;
1589                                    }
1590                                }
1591                            }//for d< N1D.size
1592                            for(Int_t d=0;d<N2D.size();d++){
1593                                Rmin = opt.ED[opt.Hyst2D[N2D[d]].HI.NED].Emin;
1594                                Rmax = opt.ED[opt.Hyst2D[N2D[d]].HI.NED].Emax;
1595                                Detrmin = opt.ED[opt.Hyst2D[N2D[d]].HI.NED].detrmin;
1596                                Detrmax = opt.ED[opt.Hyst2D[N2D[d]].HI.NED].detrmax;
1597                                Double_t Argv = 0;
1598                                Double_t LTI = 0;
1599                                Double_t DTI = 0;
1600                                if(rigev>Rmin && rigev<=Rmax && detr>Detrmin && detr<=Detrmax){
1601                                    Int_t ExpType = opt.Hyst2D[N2D[d]].HI.ExpType;
1602                                    if(ExpType == 1 || ExpType == 3) LTI = 0.16*pam_events->GetTrigLevel2()->dltime[0];
1603                                    if(ExpType == 2 || ExpType == 3) DTI = 0.01*pam_events->GetTrigLevel2()->dltime[1];
1604                                    switch  ( ExpType ){
1605                                        case 1:{Argv = LTI;break;}
1606                                        case 2:{Argv = DTI;break;}
1607                                        case 3:{Argv = LTI+DTI;break;}
1608                                        case 4:{
1609                                                Argv = opt.Hyst2D[N2D[d]].HI.LieveTime;
1610                                                opt.Hyst2D[N2D[d]].HI.LieveTime = 0;
1611                                                break;
1612                                                }
1613                                        case 5:{
1614                                                Argv = opt.Hyst2D[N2D[d]].HI.DeathTime;
1615                                                opt.Hyst2D[N2D[d]].HI.DeathTime = 0;
1616                                                break;
1617                                                }
1618                                        case 6:{
1619                                                Argv = opt.Hyst2D[N2D[d]].HI.LieveTime+opt.Hyst2D[N2D[d]].HI.DeathTime;
1620                                                opt.Hyst2D[N2D[d]].HI.LieveTime = 0;
1621                                                opt.Hyst2D[N2D[d]].HI.DeathTime = 0;
1622                                                break;
1623                                                }
1624                                    };
1625                                    Double_t Xfill = 0;
1626                                    Double_t Yfill = 0;
1627                                    //cout<<"YType = "<<opt.Hyst2D[d].YType<<endl;
1628                                    //cout<<"XType = "<<opt.Hyst2D[d].XType<<endl;
1629                                    switch ( opt.Hyst2D[N2D[d]].XType ){
1630                                        case 1:{ Xfill = TMath::RadToDeg()*PamAzim;break; }
1631                                        case 2:{ Xfill = TMath::RadToDeg()*PamZenit;break; }
1632                                        case 3:{ Xfill = SB; break; }
1633                                        case 4:{ Xfill = detr; break; }
1634                                        case 5:{ Xfill = rigev; break; }
1635                                        case 6:{ Xfill = Lat; break; }
1636                                        case 7:{ Xfill = Lon; break; }
1637                                        case 8:{ Xfill = L; break; }
1638                                        case 9:{ Xfill = Babs; break; }
1639                                        case 10:{Xfill = Alt; break; }
1640                                        case 11:{Xfill = LieveTime; break; }
1641                                        case 12:{Xfill = DeathTime; break; }
1642                                        case 13:{Xfill = SC; break; }
1643                                        case 14:{Xfill = SA; break; }
1644                                        case 17:{Xfill = log10(rigev); break;}
1645                                    };
1646                                    switch ( opt.Hyst2D[N2D[d]].YType ){
1647                                        case 1:{ Yfill = TMath::RadToDeg()*PamAzim;break; }
1648                                        case 2:{ Yfill = TMath::RadToDeg()*PamZenit;break; }
1649                                        case 3:{ Yfill = SB; break; }
1650                                        case 4:{ Yfill = detr; break; }
1651                                        case 5:{ Yfill = rigev; break; }
1652                                        case 6:{ Yfill = Lat; break; }
1653                                        case 7:{ Yfill = Lon; break; }
1654                                        case 8:{ Yfill = L; break; }
1655                                        case 9:{ Yfill = Babs; break; }
1656                                        case 10:{Yfill = Alt; break; }
1657                                        case 11:{Yfill = LieveTime; break; }
1658                                        case 12:{Yfill = DeathTime; break; }
1659                                        case 13:{Yfill = SC; break; }
1660                                        case 14:{Yfill = SA; break; }
1661                                        case 17:{Yfill = log10(rigev); break;}
1662                                    };
1663                                    if(verb){
1664                                        cout<<endl;
1665                                        cout<<"rigev = "<<rigev<<"\tdetr = "<<detr<<endl;
1666                                        cout<<"Babs = "<<Babs<<"\tL = "<<L<<endl;
1667                                        cout<<"I Fill "<<opt.Hyst2D[N2D[d]].HI.filename<<" hystogramm"<<endl;
1668                                        cout<<"Xfill = "<<Xfill<<endl<<" Yfill = "<<Yfill<<endl;
1669                                    }
1670                                    opt.Hyst2D[N2D[d]].Hyst2DF.Fill(Xfill,Yfill);
1671                                    if(ExpType!=0) opt.Hyst2D[N2D[d]].Hyst2DFExp.Fill(Xfill,Yfill,Argv);
1672                                    if(verb){
1673                                        Int_t intt;
1674                                        cin>>intt;
1675                                    }
1676                                }
1677                            }//for d<N"D.size
1678                            
1679                        }//if ntr!=-1
1680                    } // k<SS.size()
1681                    if(verb)cout<<i<<endl;
1682              } // if ntr!=-1;              } // if ntr!=-1;
               
1683          } // if GetNTrack==1;          } // if GetNTrack==1;
       
1684      } //general loop      } //general loop
   
     Double_t BinContent;  
1685            
1686      /*********************************************************/      TCanvas* Canv;
1687      //Divide canvases for angular histogramms      TCanvas* Canvdiv;
1688      /*********************************************************/      TCanvas* CanvExp;
1689            TCanvas* CanvLog;
1690      if(AnglHist){      TString Ext[3];
1691            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);  
1692            
1693      }      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  
     /**************************************************************************************/  
1694            
1695      if(PamEff){      gStyle->SetPalette(1);
       
     c7->Divide(1,2);  
       
     }  
1696            
1697      /**************************************************************************************/      for(Int_t i = 0; i<opt.Hyst1D.size(); i++){
1698      //Draw Angular Histogramms          TH1F hdiv;
1699      /**************************************************************************************/          TH1F h1n1D; TH1F h2n1D; TH1F hlog;
1700            if(opt.Hyst1D[i].HI.ExpType!=0){
1701      if(AnglHist){              hlog.SetBins(80,-1.,3.);
1702                h1n1D.SetBins(opt.Hyst1D[i].Hyst1DF.GetNbinsX(),opt.Hyst1D[i].Hyst1DF.GetXaxis()->GetXmin(),opt.Hyst1D[i].Hyst1DF.GetXaxis()->GetXmax());
1703      c1->cd(1);              h2n1D.SetBins(opt.Hyst1D[i].Hyst1DF.GetNbinsX(),opt.Hyst1D[i].Hyst1DF.GetXaxis()->GetXmin(),opt.Hyst1D[i].Hyst1DF.GetXaxis()->GetXmax());
1704      MainAxisPamelaExpositionBrazil->Draw();              //h2div.SetBins(h12D->GetNbinsX(),h12D->GetXaxis()->GetXmin(),h12D->GetXaxis()->GetXmax(),h12D->GetNbinsY(),h12D->GetYaxis()->GetXmin(),h12D->GetYaxis()->GetXmax());
1705      c1->cd(2);              for(Int_t ii = 0; ii < opt.Hyst1D[i].Hyst1DF.GetNbinsX(); ii++){
1706      MainAxisPamelaCountBrazil->Draw();                  h1n1D.SetBinContent(ii,opt.Hyst1D[i].Hyst1DF.GetBinContent(ii)/opt.Hyst1D[i].Hyst1DF.GetEntries());
1707      c1->cd(3);                  h2n1D.SetBinContent(ii,opt.Hyst1D[i].Hyst1DFExp.GetBinContent(ii)/opt.Hyst1D[i].Hyst1DFExp.GetEntries());
1708      PitchExpositionNoOrientationBrazil->Draw();              }
1709      c1->cd(4);              hdiv.SetName(opt.Hyst1D[i].Hyst1DF.GetName());
1710      PitchCountNoOrientationBrazil->Draw();              hdiv.SetBins(opt.Hyst1D[i].Hyst1DF.GetNbinsX(),opt.Hyst1D[i].Hyst1DF.GetXaxis()->GetXmin(),opt.Hyst1D[i].Hyst1DF.GetXaxis()->GetXmax());
1711      c1->cd(5);              hdiv.Divide(&h1n1D,&h2n1D);
1712      DiferenceBetweenNoOrientationAndWithOrientationBrazil->Divide(PitchExpositionNoOrientationBrazil,MainAxisPamelaExpositionBrazil);              cout<<"i = "<<i<<endl;
1713      DiferenceBetweenNoOrientationAndWithOrientationBrazil->Draw();              if(opt.Hyst1D[i].XType==17) hlog.Divide(&opt.Hyst1D[i].Hyst1DF,&binw);
1714      c1->cd(6);          }
1715      for(Int_t iu=0;iu<360;iu++){          if(opt.Hyst1D[i].HI.TypeExt > 0){
1716          BinContent = DiferenceBetweenNoOrientationAndWithOrientationBrazil->GetBinContent(iu);              Canv = new TCanvas("Canv",opt.Hyst1D[i].Hyst1DF.GetName(),1200,800);
1717          DifferenceCountInBrazil->Fill(BinContent);              Canv->cd();
1718      }              opt.Hyst1D[i].Hyst1DF.Draw();
1719      DifferenceCountInBrazil->Draw();              Canv->SaveAs((TString)outdir+"ParDirCalc_"+opt.Hyst1D[i].HI.filename+".png");
1720      //DifferenceCountInBrazil->SaveAs((TString)outdir+"PIC001A_BrazilExpositionAndCount.root");              if(opt.Hyst1D[i].HI.ExpType!=0){
1721      c1->SaveAs((TString)outdir+"PIC001A_BrazilExpositionAndCount.jpg");                  CanvExp = new TCanvas("CanvExp",opt.Hyst1D[i].Hyst1DFExp.GetName(),1200,800);
1722                        CanvExp->cd();
1723      c2->cd(1);                  opt.Hyst1D[i].Hyst1DFExp.Draw();
1724      MainAxisPamelaExpositionWithoutBrazil->Draw();                  CanvExp->SaveAs((TString)outdir+"ParDirCalc_"+opt.Hyst1D[i].HI.filename+"_Exp.png");
1725      c2->cd(2);                  Canvdiv = new TCanvas("Canvdiv",hdiv.GetName(),1200,800);
1726      MainAxisPamelaCountWithoutBrazil->Draw();                  Canvdiv->cd();
1727      c2->cd(3);                  hdiv.Draw();
1728      PitchExpositionNoOrientationWithoutBrazil->Draw();                  Canvdiv->SaveAs((TString)outdir+"ParDirCalc_"+opt.Hyst1D[i].HI.filename+"DivExp.png");
1729      c2->cd(4);              }
1730      PitchCountNoOrientationWithoutBrazil->Draw();              if(opt.Hyst1D[i].XType==17){
1731      c2->cd(5);                  CanvLog = new TCanvas("CanvLog",opt.Hyst1D[i].Hyst1DF.GetName(),1200,800);
1732      DiferenceBetweenNoOrientationAndWithOrientationWithoutBrazil->Divide(PitchExpositionNoOrientationWithoutBrazil,MainAxisPamelaExpositionWithoutBrazil);                  CanvLog->cd();
1733      DiferenceBetweenNoOrientationAndWithOrientationWithoutBrazil->Draw();                  hlog.Draw();
1734      c2->cd(6);                  CanvLog->SaveAs((TString)outdir+"ParDirCalc_"+opt.Hyst1D[i].HI.filename+"_DivBinWide.png");
1735      for(Int_t iu=0;iu<360;iu++){              }
1736          BinContent = DiferenceBetweenNoOrientationAndWithOrientationWithoutBrazil->GetBinContent(iu);          }
1737          DifferenceCountWithoutBrazil->Fill(BinContent);          if(opt.Hyst1D[i].HI.TypeExt == 0 || opt.Hyst1D[i].HI.TypeExt == 2){
1738      }              TFile fg((TString)outdir+"ParDirCalc_"+opt.Hyst1D[i].HI.filename+".root","recreate");
1739      DifferenceCountWithoutBrazil->Draw();              opt.Hyst1D[i].Hyst1DF.Write();
1740      c2->SaveAs((TString)outdir+"PIC003A_WithoutBrazilExpositionAndCount.jpg");              fg.Close();
1741                    TFile fe((TString)outdir+"ParDirCalc_"+opt.Hyst1D[i].HI.filename+"_DivWide.root","recreate");
1742      c3->cd(1);              hlog.Write();
1743      MainAxisPamelaExpositionEquator->Draw();              fe.Close();
1744      c3->cd(2);              if(opt.Hyst1D[i].HI.ExpType!=0){
1745      MainAxisPamelaCountEquator->Draw();                  TFile fgExp((TString)outdir+"ParDirCalc_"+opt.Hyst1D[i].HI.filename+"_Exp.root","recreate");
1746      c3->cd(3);                  opt.Hyst1D[i].Hyst1DFExp.Write();
1747      PitchExpositionNoOrientationEquator->Draw();                  fg.Close();
1748      c3->cd(4);                  TFile fgdiv((TString)outdir+"ParDirCalc_"+opt.Hyst1D[i].HI.filename+"_DivExp.root","recreate");
1749      PitchCountNoOrientationEquator->Draw();                  hdiv.Write();
1750      c3->cd(5);                  fgdiv.Close();
1751      DiferenceBetweenNoOrientationAndWithOrientationEquator->Divide(PitchExpositionNoOrientationEquator,MainAxisPamelaExpositionEquator);              }
1752      DiferenceBetweenNoOrientationAndWithOrientationEquator->Draw();              
1753      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");  
   
1754      }      }
1755        for(Int_t i = 0; i<opt.Hyst2D.size(); i++){
1756      /***************************************************************************************/          TH2F hdiv;
1757      //Draw hystogramms for diferences in results between calculating using DoTrack2          TH2F h1n2D; TH2F h2n2D;
1758      //and using axv[0], ayv[0]          cout<<"i = "<<i<<endl;
1759      /**************************************************************************************/          if(opt.Hyst2D[i].HI.ExpType!=0){
1760                //cout<<"Set hdiv"<<endl;
1761      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());
1762                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());
1763      c6->cd(1);              //h2div.SetBins(h12D->GetNbinsX(),h12D->GetXaxis()->GetXmin(),h12D->GetXaxis()->GetXmax(),h12D->GetNbinsY(),h12D->GetYaxis()->GetXmin(),h12D->GetYaxis()->GetXmax());
1764      DifDoTr2ANdaxvAzim->Draw();              for(Int_t ii = 0; ii < opt.Hyst2D[i].Hyst2DF.GetNbinsX(); ii++){
1765      c6->cd(2);                  for(Int_t j = 0; j < opt.Hyst2D[i].Hyst2DF.GetNbinsY(); j++){
1766      DifDoTr2ANdaxvZenit->Draw();                      h1n2D.SetBinContent(ii,j,opt.Hyst2D[i].Hyst2DF.GetBinContent(ii,j)/opt.Hyst2D[i].Hyst2DF.GetEntries());
1767      c6->SaveAs((TString)outdir+"PIC001Diference.jpg");                      h2n2D.SetBinContent(ii,j,opt.Hyst2D[i].Hyst2DFExp.GetBinContent(ii,j)/opt.Hyst2D[i].Hyst2DFExp.GetEntries());
1768                        }
1769                }
1770                hdiv.SetName(opt.Hyst2D[i].Hyst2DF.GetName());
1771                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());
1772                //cout<<"I try to divide it!!!"<<endl;
1773                hdiv.Divide(&h1n2D,&h2n2D);
1774                //cout<<"NbinsX1 = "<<h1n2D.GetNbinsX()<<" NBinsY1 = "<<h1n2D.GetNbinsY()<<endl;
1775                //cout<<"NbinsX2 = "<<h2n2D.GetNbinsX()<<" NBinsY2 = "<<h2n2D.GetNbinsY()<<endl;
1776                //cout<<"I has divide it!!! :))))\n0_o"<<endl;
1777            }
1778            //cout<<"opt.Hyst2D["<<i<<"].HI.TypeExt = "<<opt.Hyst2D[i].HI.TypeExt<<endl;
1779            if(opt.Hyst2D[i].HI.TypeExt > 0){
1780                //cout<<"opt.Hyst2D[i].HI.TypeExt = "<<opt.Hyst2D[i].HI.TypeExt<<endl;
1781                Canv = new TCanvas("Canv",opt.Hyst2D[i].Hyst2DF.GetName(),1200,800);
1782                Canv->cd();
1783                //cout<<"Canvas..."<<endl;
1784                opt.Hyst2D[i].Hyst2DF.Draw("colz");
1785                //cout<<"Draw hyst..."<<endl;
1786                Canv->SaveAs((TString)outdir+"ParDirCalc_"+opt.Hyst2D[i].HI.filename+".png");
1787                cout<<"Save canvas..."<<endl;
1788                if(opt.Hyst2D[i].HI.ExpType!=0){
1789                    CanvExp = new TCanvas("CanvExp",opt.Hyst2D[i].Hyst2DFExp.GetName(),1200,800);
1790                    CanvExp->cd();
1791                    opt.Hyst2D[i].Hyst2DFExp.Draw("colz");
1792                    CanvExp->SaveAs((TString)outdir+"ParDirCalc_"+opt.Hyst2D[i].HI.filename+"_Exp.png");
1793                    Canvdiv = new TCanvas("Canvdiv",hdiv.GetName(),1200,800);
1794                    Canvdiv->cd();
1795                    hdiv.Draw("colz");
1796                    Canvdiv->SaveAs((TString)outdir+"ParDirCalc_"+opt.Hyst2D[i].HI.filename+"DivExp.png");
1797                }
1798            }
1799            if(opt.Hyst2D[i].HI.TypeExt == 0 || opt.Hyst2D[i].HI.TypeExt == 2){
1800                TFile fg((TString)outdir+"ParDirCalc_"+opt.Hyst2D[i].HI.filename+".root","recreate");
1801                opt.Hyst2D[i].Hyst2DF.Write();
1802                fg.Close();
1803                if(opt.Hyst2D[i].HI.ExpType!=0){
1804                    TFile fgExp((TString)outdir+"ParDirCalc_"+opt.Hyst2D[i].HI.filename+"_Exp.root","recreate");
1805                    opt.Hyst2D[i].Hyst2DFExp.Write();
1806                    fg.Close();
1807                    cout<<"I'm here!"<<endl;
1808                    TFile fgdiv((TString)outdir+"ParDirCalc_"+opt.Hyst2D[i].HI.filename+"_DivExp.root","recreate");
1809                    hdiv.Write();
1810                    fgdiv.Close();
1811                    cout<<"I'm here!"<<endl;
1812                }
1813            }
1814            //delete Canv;
1815            //delete Canvdiv;
1816            //delete CanvExp;
1817            //delete CanvLog;
1818      }      }
       
     /***********************************************************************/  
     //Draw hystogram to count how much time Pamela have each Pitch-Angle  
     /***********************************************************************/  
       
     if(PamAngTime) {  
   
     c8->cd();  
     PitchTime->Draw();  
     c8->SaveAs((TString)outdir+"PIC002PamAngTime.jpg");  
1819    
     }  
       
     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");  
       
     }  
       
1820      /**************************************************************************************/      /**************************************************************************************/
1821      //Draw hystogramms for protons. Volodia      //Draw hystogramms for protons. Volodia
1822      /**************************************************************************************/      /**************************************************************************************/
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.3
  ViewVC Help
Powered by ViewVC 1.1.23