//Program selekcjonujacy przypadki typu NC i CC eksperymentu MINOS
//niezbedne pliki dodatkowe: MyClass.h oraz MyClass.C

//Autor: Krzysztof W. Fornalski, WF PW, fornalski@knf.pw.edu.pl

//wywolywanie pod srodowiskiem LOON (nakladka na ROOTa): .x moj_plik.C
//wywolywanie bezposrednie w konsoli: loon moj_plik.C

void moj_plik()
{  
  //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  //PARAMETRY WEJCIOWE:
  //////////////////////////////////////
  
  //wybierz nr zmiennej:
  int zmiennaN_1=3; //nr zmiennej do histogramow N-1
  int zmiennaPO=28; //nr zmiennej do histogramow po poprzednich cieciach
  int zmiennaTEST=97; //nr zmiennej do rozkladu testowego TEST1 i TEST2

  //ktore histogramy wlaczyc?
  bool PurEff=0; //czy wlaczyc histogram Pur i Eff (TRUE=TAK=1; FALSE=NIE=0)
  bool N_1=0; //czy wlaczyc histogram N-1 (TAK=1; NIE=0)
  bool GeV=0; //czy wlaczyc histogram rozkladu energii wiazki w GeVach (TAK=1; NIE=0)
  bool TEST1=0; //czy wlaczyc histogram testowy nr 1 (TAK=1; NIE=0) - po cieciach (=czastka)
  bool TEST2=1; //czy wlaczyc histogram testowy nr 2 (TAK=1; NIE=0) - bez ciec (=iaction)
  bool H_STP=0; //czy wlaczyc histogram z galezia STP (TAK=1; NIE=0)
  bool H_PO=0; //czy wlaczyc histogram zmiennej "zmiennaPO" po poprzednich cieciach

  //ktore zmienne usunac z selekcji? TRUE=1=zostawic; FALSE=0=usunac
  bool zmienna28=1;
  bool zmienna3=1;
  bool zmienna18=0;
  bool zmienna26=0;
  bool zmienna7=1;
  bool zmienna2=0;
  bool zmienna13=0;
  bool zmienna20=1;
  bool zmienna8=0;
  bool zmienna6=1;

  bool czyNCjakoNC=1; // do liczenia N-1; 1=(czy szukac NCjakoNC oraz CCjakoCC), czy 0=(NCjakoCC i CCjakoNC)?
  //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////


//&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
////////////////////////////////////////////////////////////////////////////////
/////////////////////////FUNKCJE WEWNETRZNE/////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
//&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

void ustawienia()  //funkcja odpowiada za ustawienia wygladu prezentowanych histogramow
{
  gStyle->SetOptStat(10000011);  // no statistics
  //  gStyle->SetOptStat(0);  //no statistics _or_
  //gStyle->SetOptStat(11111111);
  // -->all statistics (the more 1's you put,the more statistics you get)
  //  gStyle->SetOptLogy();  //log. scale !!!!!!!!!! Czyli skala logarytmiczna

  gStyle->SetOptTitle(0);  //no title
  //gStyle->SetPalette(1);   // Better color palette
  gStyle->SetCanvasColor(kWhite);     // background is no longer mouse-dropping
  gStyle->SetPalette(1,0);            // blue to red false color palette. Use 9
  //gStyle->SetFillColor(0);  //white fill color
  gStyle->SetFrameBorderMode(0);  //no frame border
  gStyle->SetCanvasBorderMode(0);  //no canvas border
  gStyle->SetPadBorderMode(0);

  // For publishing:
  gStyle->SetTextSize(0.01);
  gStyle->SetTitleSize(0.04,"xy");
  gStyle->SetTitleOffset(1.2,"x");
  gStyle->SetTitleOffset(1.2,"y");
  //testowa=1;
}
//////////////////////////////////////////////////////////////////////////////////

float sqr(float xxx) //kwadrat
{
  return xxx*xxx;
}
//////////////////////////////////////////////////////////////////////////////////

void hist_TEST() //Wyswietla histogram testujacy rozklady zmiennych dla NC/CC (po cieciach)
{
  //tworzenie okna histogramu
  //THStack hs ("hs","test stack;Badana zmienna;Liczba przypadkow");
  TCanvas* TEST = new TCanvas ("canvas-testowy","TEST",70,300,1200,600);

  //stara wersja - dwa okna
  /*
  // TCanvas *c0 = new TCanvas("eff", "Efektywnosc", 0, 10, 800, 600);
  //TEST->Divide(2,1);
  //TEST->cd(1);
  TEST_NC->Draw();
  //TEST->cd(2);
  TEST_CC->Draw("SAME");
  c0->Update();
  c0->Modified();*/

  //TCanvas* TEST = new TCanvas ("canvas-testowy","Test zmiennej NC/CC",70,300,1200,600);
  //TEST_CC->SetFillColor(kBlue); //niebieski = CC
  //TEST_NC->SetFillColor(kGreen); //zielony =NC
  //hs.Add(TEST_NC);
  //hs.Add(TEST_CC);
  //TEST-CC->SetPalette(kBlue);
  //TEST-NC->SetPalette(kGreen);

  //nakladanie na siebie w 2 kolorach w JEDNNYM OKNIE - to co na seminarium pokazywalem:

  TEST_NC->SetFillColor (3);
  TEST_NC->SetFillStyle (3001);
  TEST_CC->SetFillColor (4);
  TEST_CC->SetFillStyle (3001);

  gStyle->SetOptStat(0);
  gStyle->SetTitleOffset(1.2,"x");
  gStyle->SetTitleOffset(1.2,"y");

  //TEST->GetXaxis()->SetTitle("Badana zmienna");
  //TEST->GetYaxis()->SetTitle("Ilosc przypadkow");
  TEST->Divide (1,1);
  TEST->cd(1);
  //hs.Draw();
  TEST_CC->Draw();  //w przypadku wyjscia poza skale, trzeba zamienic NC z CC
  TEST_NC->Draw("SAME");
  TEST->Update();
  TEST->Modified();

  TLegend *legend = new TLegend (0.8, 0.88, 0.88, 0.7);
  legend->AddEntry (TEST_NC, "NC", "f");
  legend->AddEntry (TEST_CC, "CC", "f");
  legend->Draw ();
}
//////////////////////////////////////////////////////////////////////////////////

void hist_TEST2() //dzialanie analogiczne jak hist_TEST, lecz przed cieciami. Druga wersja potrzebna byla dla porownania niektorych danych (po przejsciu N ciec oraz N-1)
{
  //THStack hs2 ("hs2","test stack2;Badana zmienna;Liczba przypadkow");
  TCanvas* TEST2 = new TCanvas ("canvas-testowy2","TEST 2",70,300,1200,600);

  TEST_NC2->SetFillColor (3);
  TEST_NC2->SetFillStyle (3001);
  TEST_CC2->SetFillColor (4);
  TEST_CC2->SetFillStyle (3001);

  gStyle->SetOptStat(0);
  gStyle->SetTitleOffset(1.2,"x");
  gStyle->SetTitleOffset(1.2,"y");

  //TEST->GetXaxis()->SetTitle("Badana zmienna");
  //TEST->GetYaxis()->SetTitle("Ilosc przypadkow");
  TEST2->Divide (1,1);
  TEST2->cd(1);
  //hs.Draw();
  TEST_CC2->Draw();
  TEST_NC2->Draw("SAME");
  TEST2->Update();
  TEST2->Modified();

  TLegend *legend4 = new TLegend (0.8, 0.88, 0.88, 0.7);
  legend4->AddEntry (TEST_NC2, "NC", "f");
  legend4->AddEntry (TEST_CC2, "CC", "f");
  legend4->Draw ();
}
//////////////////////////////////////////////////////////////////////////////////

void HIST_PO (TH1D* PO_NC, TH1D* PO_CC)
{
  TCanvas* PO = new TCanvas ("PO","PO",70,300,1200,600);

  PO_NC->SetFillColor (3);
  PO_NC->SetFillStyle (3001);
  PO_CC->SetFillColor (4);
  PO_CC->SetFillStyle (3001);

  gStyle->SetOptStat(0);
  gStyle->SetTitleOffset(1.2,"x");
  gStyle->SetTitleOffset(1.2,"y");

  PO->Divide (1,1);
  PO->cd(1);
  PO_CC->Draw();
  PO_NC->Draw("SAME");
  PO->Update();
  PO->Modified();

  TLegend *legend4 = new TLegend (0.8, 0.88, 0.88, 0.7);
  legend4->AddEntry (PO_NC, "NC", "f");
  legend4->AddEntry (PO_CC, "CC", "f");
  legend4->Draw ();
}
//////////////////////////////////////////////////////////////////////////////////

void Nmniej1(TH1F* N1_wszystkoNC, TH1F* N1_1ciecieNC, TH1F* N1_cieciaNC, TH1F* N1_wszystkoCC, TH1F* N1_1ciecieCC, TH1F* N1_cieciaCC) //Wyswietla oba histogramy (NC i CC) obrazujacy liczbe przypadkow po N-1 ciec 
{
  //NC:
  //THStack hsN1 ("hs2","test stack 2;Badana zmienna;Liczba przypadkow");
  TCanvas* Nmniej1NC = new TCanvas ("canvas-N-1 NC","N-1 NC",70,300,1200,600);
  N1_wszystkoNC->SetLineColor (2);
  N1_wszystkoNC->SetLineWidth (3);
//  N1_wszystkoNC->SetFillStyle (4000);
  N1_1ciecieNC->SetLineColor (8);
  N1_1ciecieNC->SetLineStyle (4);
  N1_1ciecieNC->SetLineWidth (3);
//  N1_1ciecieNC->SetFillStyle (3001);
  N1_cieciaNC->SetLineColor (9);
  N1_cieciaNC->SetLineStyle (2);
  N1_cieciaNC->SetLineWidth (3);
//  N1_cieciaNC->SetFillStyle (3001);

  gStyle->SetOptStat(0); //Opt
  gStyle->SetTitleOffset(1.2,"x");
  gStyle->SetTitleOffset(1.2,"y");

  Nmniej1NC->Divide (1,1);
  Nmniej1NC->cd(1);
  //hs.Draw();
  N1_wszystkoNC->Draw();
  N1_cieciaNC->Draw("SAME");
  N1_1ciecieNC->Draw("SAME");
  Nmniej1NC->Update();
  Nmniej1NC->Modified();

/*  N1_wszystkoNC->SetTitle("");
  N1_wszystkoNC->SetXTitle("Energia [GeV]");
  N1_wszystkoNC->SetYTitle("Eff - efektywnosc NC");
  N1_wszystkoNC->SetStats(0);*/

  TLegend *legend2 = new TLegend (0.8, 0.88, 0.88, 0.7);
  legend2->AddEntry (N1_wszystkoNC, "pierwotna wiazka", "f");
  legend2->AddEntry (N1_1ciecieNC, "N-1 ciec", "f");
  legend2->AddEntry (N1_cieciaNC, "wszystkie ciecia", "f");
  legend2->Draw ();
  
  //CC:
  //THStack hsN2 ("hs3","test stack 3;Badana zmienna;Liczba przypadkow");
  TCanvas* Nmniej1CC = new TCanvas ("canvas-N-1 CC","N-1 CC",70,300,1200,600);
  N1_wszystkoCC->SetLineColor (2);
  N1_wszystkoCC->SetLineWidth (3);
//  N1_wszystkoNC->SetFillStyle (4000);
  N1_1ciecieCC->SetLineColor (8);
  N1_1ciecieCC->SetLineStyle (4);
  N1_1ciecieCC->SetLineWidth (3);
//  N1_1ciecieNC->SetFillStyle (3001);
  N1_cieciaCC->SetLineColor (9);
  N1_cieciaCC->SetLineStyle (2);
  N1_cieciaCC->SetLineWidth (3);

  gStyle->SetOptStat(0);
  gStyle->SetTitleOffset(1.2,"x");
  gStyle->SetTitleOffset(1.2,"y");

  Nmniej1CC->Divide (1,1);
  Nmniej1CC->cd(1);
  //hs.Draw();

  N1_wszystkoCC->Draw();
  N1_cieciaCC->Draw("SAME");
  N1_1ciecieCC->Draw("SAME");
  Nmniej1CC->Update();
  Nmniej1CC->Modified();

  TLegend *legend3 = new TLegend (0.8, 0.88, 0.88, 0.7);
  legend3->AddEntry (N1_wszystkoCC, "pierwotna wiazka", "f");
  legend3->AddEntry (N1_1ciecieCC, "N-1 ciec", "f");
  legend3->AddEntry (N1_cieciaCC, "wszystkie ciecia", "f");
  legend3->Draw ();
}
//////////////////////////////////////////////////////////////////////////////////

void pelne_energie() //rysuje 4 histogramy energetyczne w jednym oknie
{
  //tylko tutaj LOGY
  TCanvas* canvas = new TCanvas ("canvas","Histogramy- pelne energie",30,20,1200,900);
  canvas->Divide(2,2);
  canvas->cd(1);
  hist_ph_pe->Draw();
  canvas->cd(2);
  hist_ph_raw->Draw();
  canvas->cd(3);
  hist_ph_siglin->Draw();
  canvas->cd(4);
  hist_ph_sigcor->Draw();
  canvas->Update();
  canvas->Modified();
}
//////////////////////////////////////////////////////////////////////////////////

void niskie_energie() //rysuje 2 histogramy dla niskich energii
{
  TCanvas* canvasLOW_CC = new TCanvas ("canvas2","Histogramy- niskie energie CC",50,50,1200,900);
  TCanvas* canvasLOW_NC = new TCanvas ("canvas2a","Histogramy- niskie energie NC",50,50,1200,900);

  canvasLOW_CC->Divide(2,2);
  canvasLOW_CC->cd(1);
  hist_ph_peLOW_CC->Draw();
  canvasLOW_CC->cd(2);
  hist_ph_rawLOW_CC->Draw();
  canvasLOW_CC->cd(3);
  hist_ph_siglinLOW_CC->Draw();
  canvasLOW_CC->cd(4);
  hist_ph_sigcorLOW_CC->Draw();
  canvasLOW_CC->Update();
  canvasLOW_CC->Modified();

  canvasLOW_NC->Divide(2,2);
  canvasLOW_NC->cd(1);
  hist_ph_peLOW_NC->Draw();
  canvasLOW_NC->cd(2);
  hist_ph_rawLOW_NC->Draw();
  canvasLOW_NC->cd(3);
  hist_ph_siglinLOW_NC->Draw();
  canvasLOW_NC->cd(4);
  hist_ph_sigcorLOW_NC->Draw();
  canvasLOW_NC->Update();
  canvasLOW_NC->Modified();
}
//////////////////////////////////////////////////////////////////////////////////

void hist_NC_CC() //histogram rozdzielajacy przypadki NC od CC dla mc.iaction
{
  //glowny rozdzial przypadkow NC i CC
  TCanvas* NC_CC = new TCanvas ("canvas3","Rozdzial przypadkow NC i CC",70,60,1200,900);
  NC_CC->Divide(2,2);
  NC_CC->cd(1);
  sigcor_CC->Draw();
  NC_CC->cd(2);
  sigcor_NC->Draw();
  NC_CC->cd(3);
  sigcor_UNKNOWN->Draw();
  NC_CC->cd(4);
  sigcor_Suma->Draw();
  NC_CC->Update();
  NC_CC->Modified();
}
//////////////////////////////////////////////////////////////////////////////////

void hist_ZERO_JEDEN(TH1D* sigcor_JEDEN, TH1D* sigcor_ZERO) //rysuje histogram z rozkladem energii w zaleznosci od NC/CC (mc.iaction)
{
  TCanvas* ZERO_JEDEN = new TCanvas ("canvas4","mc_iaction",70,270,1200,600);
  ZERO_JEDEN->Divide(2,1);
  ZERO_JEDEN->cd(1);
  sigcor_JEDEN->Draw();
  ZERO_JEDEN->cd(2);
  sigcor_ZERO->Draw();
  ZERO_JEDEN->Update();
  ZERO_JEDEN->Modified();
}
//////////////////////////////////////////////////////////////////////////////////

void Pur_Eff(TH1F* CCjakoCC, TH1F* CCjakoNC, TH1F* NCjakoNC, TH1F* NCjakoCC, TH1F* CCall, TH1F* NCall, float itCCjakoCC, float itCCjakoNC, float itNCjakoNC, float itNCjakoCC, float iIleCC, float iIleNC) //tworzy 4 histogramy odpowiadajace za czystosc (Pur) i efektywnosc (Eff) dla NC i CC
{
  TH1F *hPurCC = new TH1F("PurCC","PurCC",50,0,10);
  TH1F *hEffCC = new TH1F("EffCC","EffCC",50,0,10);    
  TH1F *hPurNC = new TH1F("PurNC","PurNC",50,0,10);
  TH1F *hEffNC = new TH1F("EffNC","EffNC",50,0,10);

  //Obliczenia
  TH1F *temp = new TH1F("temp","temp;Energia [GeV];CC_jako_CC + NC_jako_CC",50,0,10); //powyrzucac opisy osi
  temp -> Add(CCjakoCC,NCjakoCC,1,1);
  hPurCC -> Divide(CCjakoCC,temp,1,1);

  hEffCC -> Divide(CCjakoCC, CCall, 1, 1);
      
  TH1F *temp1 = new TH1F("temp1","temp1;Energia [GeV];CC_jako_NC + NC_jako_NC",50,0,10);
  temp1 -> Add(NCjakoNC,CCjakoNC,1,1);
  hPurNC -> Divide(NCjakoNC,temp1,1,1);
  
  hEffNC -> Divide(NCjakoNC, NCall, 1, 1);
  //histogramy zmiennych dzielacych (wewnetrznych)
  /*TCanvas *c0 = new TCanvas("effNC", "Eff - efektywnosc NC", 0, 10, 800, 600);
  temp1->Draw();
  c0->Update();
  c0->Modified();*/

  ///////////Eff NC:
  TCanvas *c0 = new TCanvas("effNC", "Eff - efektywnosc NC", 0, 10, 800, 600);
  hEffNC->SetMinimum(0);
  hEffNC->SetMaximum(1);

  //bledy:
/*  Int_t nbinsx = hEffNC->GetNbinsX();
  for (Int_t binx=0; binx<=nbinsx+1; binx++)
    hEffNC->SetBinError(binx,sqrt(itNCjakoNC*(iIleNC-itNCjakoNC)/iIleNC) );
*/
//sqrt(sqr(sqrt(itNCjakoNC)/iIleNC)+sqr(sqrt(iIleNC)*(itNCjakoNC)/sqr(iIleNC))));
 

  //hEffCC->Draw(); //mozna dla CC i NC na jednym wykresie, wowczas przy dolnym ->Draw("SAME");
  hEffNC->Draw();  
  c0->Update();
  c0->Modified();
  hEffNC->SetTitle("");
  hEffNC->SetXTitle("Energia [GeV]");
  hEffNC->SetYTitle("Eff - efektywnosc NC");
  //hEffCC->SetStats(0);
  hEffNC->SetStats(0);
  //TLegend *legend = new TLegend (0.65, 0.72, 0.86, 0.89);
  //TLegend *legend = new TLegend (0.8, 0.14, 0.88, 0.31);
  //legend->AddEntry (hEffNC, "NC", "P");
  //legend->AddEntry (hEffCC, "CC", "P");
  //legend->Draw ();
  
  ///////////Eff CC:
  TCanvas *c0CC = new TCanvas("effCC", "Eff - efektywnosc CC", 0, 10, 800, 600);
  hEffCC->SetMinimum(0);
  hEffCC->SetMaximum(1);

  //bledy:
/*  Int_t nbinsx = hEffCC->GetNbinsX();
  for (Int_t binx=0; binx<=nbinsx+1; binx++)
    hEffCC->SetBinError(binx,sqrt(itCCjakoCC*(iIleCC-itCCjakoCC)/iIleCC) );
*/
//sqrt(sqr(sqrt(itCCjakoCC)/iIleCC)+sqr(sqrt(iIleCC)*(itCCjakoCC)/sqr(iIleCC))));


  hEffCC->Draw();
  //hEffNC->Draw("SAME");
  c0CC->Update();
  c0CC->Modified();
  hEffCC->SetTitle("");
  hEffCC->SetXTitle("Energia [GeV]");
  hEffCC->SetYTitle("Eff - efektywnosc CC");
  hEffCC->SetStats(0);
  //hEffNC->SetStats(0);
  //TLegend *legend = new TLegend (0.65, 0.72, 0.86, 0.89);
  //TLegend *legend = new TLegend (0.8, 0.14, 0.88, 0.31);
  //legend->AddEntry (hEffNC, "NC", "P");
  //legend->AddEntry (hEffCC, "CC", "P");
  //legend->Draw ();

  /////////////Pur CC:
  TCanvas *c1 = new TCanvas("purCC", "Pur - czystosc CC", 100, 110, 800, 600);    
  hPurCC->SetMinimum(0);
  hPurCC->SetMaximum(1); 

  //bledy
/*  nbinsx = hPurCC->GetNbinsX();
  for (Int_t binx=0; binx<=nbinsx+1; binx++)
    hPurCC->SetBinError(binx,sqrt(itCCjakoCC*itNCjakoCC/(itCCjakoCC+itNCjakoCC)) );
*/
//sqrt(sqr(sqrt(itCCjakoCC)*itNCjakoCC/sqr(itCCjakoCC+itNCjakoCC))+sqr(sqrt(itNCjakoCC)*(itCCjakoCC)/sqr(itCCjakoCC+itNCjakoCC))));


  hPurCC->Draw();
  //hPurNC->Draw("SAME");
  c1->Update();
  c1->Modified();
  hPurCC->SetTitle("");
  hPurCC->SetXTitle("Energia [GeV]");
  hPurCC->SetYTitle("Pur - czystosc CC");
  hPurCC->SetStats(0);
  //hPurNC->SetStats(0);
  //TLegend *legend1 = new TLegend (0.65, 0.72, 0.86, 0.89);
  //TLegend *legend1 = new TLegend (0.8, 0.14, 0.88, 0.31);
  //legend1->AddEntry (hPurNC, "NC", "P");
  //legend1->AddEntry (hPurCC, "CC", "P");
  //legend1->Draw ();

  /////////////Pur NC:
  TCanvas *c1NC = new TCanvas("purNC", "Pur - czystosc NC", 100, 110, 800, 600);
  hPurNC->SetMinimum(0);
  hPurNC->SetMaximum(1);

  //bledy
/*  nbinsx = hPurNC->GetNbinsX();
  for (Int_t binx=0; binx<=nbinsx+1; binx++)
    hPurNC->SetBinError(binx,sqrt(itNCjakoNC*itCCjakoNC/(itNCjakoNC+itCCjakoNC)) );
*/
//sqrt(sqr(sqrt(itNCjakoNC)*itCCjakoNC/sqr(itNCjakoNC+itCCjakoNC))+sqr(sqrt(itCCjakoNC)*(itNCjakoNC)/sqr(itNCjakoNC+itCCjakoNC))));

  
  hPurNC->Draw();
  //hPurNC->Draw("SAME");
  c1NC->Update();
  c1NC->Modified();
  hPurNC->SetTitle("");
  hPurNC->SetXTitle("Energia [GeV]");
  hPurNC->SetYTitle("Pur - czystosc NC");
  hPurNC->SetStats(0);
  //hPurNC->SetErrorX(3);

  hPurCC->SetMarkerStyle(22);
  hPurNC->SetMarkerStyle(23);
  hEffCC->SetMarkerStyle(22);
  hEffNC->SetMarkerStyle(23);
  hPurCC->SetMarkerColor(kBlue);
  hPurNC->SetMarkerColor(kRed);
  hEffCC->SetMarkerColor(kBlue);
  hEffNC->SetMarkerColor(kRed);
//  gStyle->SetErrorY(3);
  //c0->SaveAs("autoeff.eps"); //zapisuje histogramy do pliku
  //c1->SaveAs("autopur.eps");

/*cout<<endl<<" blad EffNC: "<<sqrt(sqr(sqrt(itNCjakoNC)/iIleNC)+sqr(sqrt(iIleNC)*(itNCjakoNC)/sqr(iIleNC)));
cout<<endl<<" blad EffCC: "<<sqrt(sqr(sqrt(itCCjakoCC)/iIleCC)+sqr(sqrt(iIleCC)*(itCCjakoCC)/sqr(iIleCC)));
cout<<endl<<" blad PurCC: "<<sqrt(sqr(sqrt(itCCjakoCC)*itNCjakoCC/sqr(itCCjakoCC+itNCjakoCC))+sqr(sqrt(itNCjakoCC)*(itCCjakoCC)/sqr(itCCjakoCC+itNCjakoCC)));
cout<<endl<<" blad PurNC: "<<sqrt(sqr(sqrt(itNCjakoNC)*itCCjakoNC/sqr(itNCjakoNC+itCCjakoNC))+sqr(sqrt(itCCjakoNC)*(itNCjakoNC)/sqr(itNCjakoNC+itCCjakoNC)))<<endl;
 */

}
//////////////////////////////////////////////////////////////////////////////////

void hist_STP(TH1D* stp_PH) //dane z galezi STP - histogram
{
  TCanvas* STP = new TCanvas ("canvas6","STP.PH0.SIGCOR & STP.PH1.SIGCOR",70,300,600,600);
  //wyniki petli po STP.*
  STP->Divide(1,1);
  STP->cd(1);
  stp_PH->Draw();
  //STP->cd(2);
  //stp_PH1->Draw();
  STP->Update();
  STP->Modified();
}
//////////////////////////////////////////////////////////////////////////////////

void galaz_MC() //2 histogramy dla CC i NC - dane z galezi MC
{
  TCanvas* MC_NC = new TCanvas ("canvas7","dane_z_galezi_MC_NC",70,320,1200,600);
  TCanvas* MC_CC = new TCanvas ("canvas7a","dane_z_galezi_MC_CC",70,320,1200,600);
  MC_NC->Divide(2,1);
  MC_NC->cd(1);
  mc_P4NEU_NC->Draw();
  MC_NC->cd(2);
  mc_SUMA_NC->Draw();
  MC_NC->Update();
  MC_NC->Modified();
 
  MC_CC->Divide(2,1);
  MC_CC->cd(1);
  mc_P4NEU_CC->Draw();
  MC_CC->cd(2);
  mc_SUMA_CC->Draw();
  MC_CC->Update();
  MC_CC->Modified();
}
//////////////////////////////////////////////////////////////////////////////////

void profile_histogram() //histogram z "best fit" - dzieki temu znajdujemy przelicznik miedzy sigcor a GeV
{
  TCanvas* Profil = new TCanvas ("canvas8","Profile histogram (dla NC)",70,320,600,600);
  //Double_t fitFunction (Double_t *x) {return x;}
  //TF1 *fitFcn = new TF1("fitFcn","x","","");
  Profil->Divide(1,1);
  Profil->cd(1);
  profil->Fit("pol1");//fitFcn","","",-2,2);
  profil->Draw();
  Profil->Update();
  Profil->Modified();
}
//////////////////////////////////////////////////////////////////////////////////

void hist_dlugosc() //drugosc przypadku/toru
{
  TCanvas* Dlugosc = new TCanvas ("canvas9","Dlugosc przypadku/toru",70,320,600,600);
  Dlugosc->Divide(1,1);
  Dlugosc->cd(1);
  dlugosc->Draw(); //"dlugosc" a nie "Dlugosc"!
  Dlugosc->Update();
  Dlugosc->Modified();
}

//////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////


//&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
////////////////////////////////////////////////////////////////////////////////
//////////////////////////PROGRAM  WLASCIWY/////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
//&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

  //gROOT->Reset();
  gSystem->Load("libCandNtupleSR.so");
  gSystem->Load("libTruthHelperNtuple.so");
  gSystem->Load("libMCNtuple.so");
  gSystem->Load("libStandardNtuple.so");
  //gSystem->Load("libEvent.so");

  gROOT->LoadMacro("MyClass.C");

  MyClass t; //podstawowa klasa

  //ustawienia wyswietlania histogramow
  gROOT->SetStyle("Plain");     
  gROOT->ForceStyle();          
  ustawienia(); //dokladne ustawienia histogramow


  //////////////////////////////////////
  //DEKLARACJE ZMIENNYCH
  //////////////////////////////////////

  Long64_t nentries = t.fChain->GetEntries(); //pobiera calkowita liczbe przypadkow z wczytanej paczki danych. Inaczej: liczba iteracji w petli glownej programu

  //CC i NC - potrzebne do Pur i Eff
  int itCCjakoCC = 0;
  int itCCjakoNC = 0;
  int itNCjakoNC = 0;
  int itNCjakoCC = 0;
  int iIleCC=0,iIleNC=0;

  //zmienne pomocnicze, przeliczniki, liczniki, etc.
  Int_t nbytes = 0, nb = 0;
  int zero=0, jeden=0;
  int licznik1=0, licznik2=0, licznik3=0, licznik_wewCC=0, licznik_wewNC=0, licznik_wewUNKNOWN=0;
  int ileN_1_CC=0, ileN_1_NC=0, ileN_NC=0, ileN_CC=0, ileALL_NC=0, ileALL_CC=0;

  //zmienne do liczenia ile przypadkow NC i CC przybylo po kazdym cieciu:
  int ileNCwszystkich=0, ileCCwszystkich=0, ileNCpo28=0, ileCCpo28=0, ileNCpo3=0, ileCCpo3=0, ileNCpo18=0, ileCCpo18=0, ileCCpo26=0, ileCCpo7bis=0, ileNCpo2=0, ileCCpo13=0, ileCCpo20=0, ileNCpo20=0, ileCCpo8=0, ileNCpo8=0, ileNCpo6=0, ileCCelse=0;


  //skalowanie wykresow:

//zmiennaPO=7;

  double od=0; //minimalna wartosc na osi X
  double skalaN=1000; //N-1 - maksymalna wartosc na osi X
  double skala=1000; //TEST
  double skala_2=1000; //PO

  const char* napisN_1="N-1;badana zmienna;liczba przypadkow";
  const char* napisTEST="TEST_X;badana zmienna;liczba przypadkow";
  const char* napisPO="PO;badana zmienna;liczba przypadkow";

  switch (zmiennaN_1)  //zakresy na osi X histogramow N-1 dla roznych zmiennych
  {
    case 28: skalaN=1.5; napisN_1="N-1;sygnal z pierwszych 25% plaszczyzn do sygnalu wiodacej kaskady;liczba przypadkow"; break;
    case 3: skalaN=50; napisN_1="N-1;liczba zapalonych plaszczyzn rejestrujacych wiecej niz 2 fotoelektrony;liczba przypadkow"; break;
    case 18: skalaN=500; napisN_1="N-1;sygnal z kaskady o najwiekszej liczbie impulsow;liczba przypadkow"; break;
    case 26: skalaN=1; napisN_1="N-1;sygnal z pierwszych 10% plaszczyzn do sygnalu wiodacej kaskady;liczba przypadkow"; break;
    case 7: skalaN=10000; napisN_1="N-1;energia przypadku sigcor do ilosci wszystkich zapalonych plaszczyzn;liczba przypadkow"; break;
    case 2: skalaN=1000; napisN_1="N-1;sygnal przypadku mierzony w fotoelektronach;liczba przypadkow"; break;
    case 13: skalaN=25; napisN_1="N-1;dlugosc wiodacego toru minus dlugosc wiodacej kaskady;liczba przypadkow"; break;
    case 20: skalaN=100000; napisN_1="N-1;sygnal wiodacej kaskady mierzony w fotoelektronach;liczba przypadkow"; break;
    case 8: skalaN=0.05; napisN_1="N-1;zasieg wiodacego toru dzielony przez energie przypadku sigcor;liczba przypadkow"; break;
    case 6: skalaN=500; napisN_1="N-1;liczba zapalonych paskow scyntylacyjnych w przypadku;liczba przypadkow"; break;
  }

  switch (zmiennaTEST)  //zakresy na osi X histogramow testowych (1 i 2) dla roznych zmiennych
  {
    case 28: skala=1.5; napisTEST="TEST;sygnal z pierwszych 25% plaszczyzn do sygnalu wiodacej kaskady;liczba przypadkow"; break;
    case 3: skala=100; napisTEST="TEST;liczba zapalonych plaszczyzn rejestrujacych wiecej niz 2 fotoelektrony;liczba przypadkow"; break;
    case 18: skala=500; napisTEST="TEST;sygnal z kaskady o najwiekszej liczbie impulsow;liczba przypadkow"; break;
    case 26: skala=1; napisTEST="TEST;sygnal z pierwszych 10% plaszczyzn do sygnalu wiodacej kaskady;liczba przypadkow"; break;
    case 7: skala=10000; napisTEST="TEST;energia przypadku sigcor do ilosci wszystkich zapalonych plaszczyzn;liczba przypadkow"; break;
    case 2: skala=5000; napisTEST="TEST;sygnal przypadku mierzony w fotoelektronach;liczba przypadkow"; break;
    case 13: skala=50; napisTEST="TEST;dlugosc wiodacego toru minus dlugosc wiodacej kaskady;liczba przypadkow"; break;
    case 20: skala=150000; napisTEST="TEST;sygnal wiodacej kaskady mierzony w fotoelektronach;liczba przypadkow"; break;
    case 8: skala=0.1; napisTEST="TEST;zasieg wiodacego toru dzielony przez energie przypadku sigcor;liczba przypadkow"; break;
    case 6: skala=1000; napisTEST="TEST;liczba zapalonych paskow scyntylacyjnych w przypadku;liczba przypadkow"; break;
    case 0: skala=10; napisTEST="TEST;Energia [GeV];liczba przypadkow"; break;
    case 91: skala=250000; napisTEST="TEST;evthdr.ph.sigcor;liczba przypadkow"; break;
    case 92: skala=2500; napisTEST="TEST;evthdr.ph.pe;liczba przypadkow"; break;
    case 93: skala=250; napisTEST="TEST;evthdr.nstrip;liczba przypadkow"; break;
    case 94: skala=250; napisTEST="TEST;evthdr.planeall.n;liczba przypadkow"; break;
    case 95: skala=250; napisTEST="TEST;evthdr.plane.n;liczba przypadkow"; break;
    case 96: skala=5; napisTEST="TEST;evthdr.ntrack;liczba przypadkow"; break;
    case 97: skala=5; napisTEST="TEST;evthdr.nshower;liczba przypadkow"; break;
  }

  switch (zmiennaPO)  //zakresy na osi X histogramow po poprzednich cieciach
  {
    case 28: skala_2=1.5; napisPO="PO;sygnal z pierwszych 25% plaszczyzn do sygnalu wiodacej kaskady;liczba przypadkow"; break;
    case 3: skala_2=50; napisPO="PO;liczba zapalonych plaszczyzn rejestrujacych wiecej niz 2 fotoelektrony;liczba przypadkow"; break;
    case 18: skala_2=500; napisPO="PO;sygnal z kaskady o najwiekszej liczbie impulsow;liczba przypadkow"; break;
    case 26: skala_2=0.5; napisPO="PO;sygnal z pierwszych 10% plaszczyzn do sygnalu wiodacej kaskady;liczba przypadkow"; break;
    case 7: skala_2=10000; napisPO="PO;energia przypadku sigcor do ilosci wszystkich zapalonych plaszczyzn;liczba przypadkow"; break;
    case 2: skala_2=2500; napisPO="PO;sygnal przypadku mierzony w fotoelektronach;liczba przypadkow"; break;
    case 13: skala_2=25; napisN_1="N-1;dlugosc wiodacego toru minus dlugosc wiodacej kaskady;liczba przypadkow"; break;
    case 20: skala_2=150000; napisPO="PO;sygnal wiodacej kaskady mierzony w fotoelektronach;liczba przypadkow"; break;
    case 8: skala_2=0.05; napisPO="PO;zasieg wiodacego toru dzielony przez energie przypadku sigcor;liczba przypadkow"; break;
    case 6: skala_2=250; napisPO="PO;liczba zapalonych paskow scyntylacyjnych w przypadku;liczba przypadkow"; break;
    case 0: skala_2=10; napisPO="PO;Przypadki niesklasyfikowane [GeV];liczba przypadkow"; break;
  }

  double CC=1, NC=0;
  if (czyNCjakoNC!=1) {CC=0; NC=1;}
  

  //////////////////////////////Zmienne z plikow:

  t.fChain.SetBranchStatus("*",0); //wylacza wszystkie zmienne (konieczne, inaczej pamiec sie zapcha)

  //aktywacja TYLKO interesujacych nas zmiennych
  t.fChain.SetBranchStatus("detsim.snarlDigits",1);
  t.fChain.SetBranchStatus("mc*",1); //cala galaz MC
  t.fChain.SetBranchStatus("evthdr*",1);
  t.fChain.SetBranchStatus("stp*",1);
  t.fChain.SetBranchStatus("mc.iaction",1);
  t.fChain.SetBranchStatus("mc.p4neu*",1);
  t.fChain.SetBranchStatus("mc.p4shw*",1);
  t.fChain.SetBranchStatus("clu.id",1);
  t.fChain.SetBranchStatus("trk.ph.sigcor",1);
  t.fChain.SetBranchStatus("trk.momentum.qp",1);
  t.fChain.SetBranchStatus("trk.range",1);
  t.fChain.SetBranchStatus("shw.plane.endu",1);
  t.fChain.SetBranchStatus("shw.plane.begu",1);
  t.fChain.SetBranchStatus("trk.plane.endu",1);
  t.fChain.SetBranchStatus("trk.plane.begu",1);     
  t.fChain.SetBranchStatus("shw.plane.endv",1);
  t.fChain.SetBranchStatus("shw.plane.begv",1);
  t.fChain.SetBranchStatus("trk.plane.endv",1);
  t.fChain.SetBranchStatus("trk.plane.begv",1);         
  t.fChain.SetBranchStatus("trk.ph.sigcor",1);
  t.fChain.SetBranchStatus("trk.plane.n",1);
  t.fChain.SetBranchStatus("trk.range",1);
  t.fChain.SetBranchStatus("trk.ds",1);
  t.fChain.SetBranchStatus("stp.plane",1);
  t.fChain.SetBranchStatus("stp.ph0.sigcor",1);
  t.fChain.SetBranchStatus("stp.ph1.sigcor",1);
  t.fChain.SetBranchStatus("shw.ph.sigcor",1);
  t.fChain.SetBranchStatus("shw.ndigit",1);
  t.fChain.SetBranchStatus("shw.ph*",1);
  t.fChain.SetBranchStatus("trk.ph*",1);
  t.fChain.SetBranchStatus("trk.fit*",1);
  t.fChain.SetBranchStatus("trk.vtx.dcosz",1);
  t.fChain.SetBranchStatus("shw.vtx.dcosz",1);
  t.fChain.SetBranchStatus("trk.plane.n",1);
  t.fChain.SetBranchStatus("trk.nstrip",1);
  t.fChain.SetBranchStatus("evt.stp",1);

  ////////////////////////Zmienne niezbedne do rysowania histogramow


  //hist_TEST:
  TH1D* TEST_NC = new TH1D ("TEST_NC",napisTEST,50.,od,skala); //3. liczba od konca=liczba binow
  TH1D* TEST_CC = new TH1D ("TEST_CC",napisTEST,50.,od,skala);

  //hist_TEST2:
  TH1D* TEST_NC2 = new TH1D ("TEST_NC2",napisTEST,50.,od,skala);
  TH1D* TEST_CC2 = new TH1D ("TEST_CC2",napisTEST,50.,od,skala);

  //Nmniej1:
  TH1F* N1_wszystkoNC = new TH1F ("N-1a NC",napisN_1,50.,od,skalaN);
  TH1F* N1_1ciecieNC = new TH1F ("N-1b NC",napisN_1,50.,od,skalaN);
  TH1F* N1_cieciaNC = new TH1F ("N-1c NC",napisN_1,50.,od,skalaN);
  TH1F* N1_wszystkoCC = new TH1F ("N-1a CC",napisN_1,50.,od,skalaN);
  TH1F* N1_1ciecieCC = new TH1F ("N-1b CC",napisN_1,50.,od,skalaN);
  TH1F* N1_cieciaCC = new TH1F ("N-1c CC",napisN_1,50.,od,skalaN);

  //PO poprzednich ciechach
  TH1D* PO_NC = new TH1D ("PO NC",napisPO,50.,od,skala_2);
  TH1D* PO_CC = new TH1D ("PO CC",napisPO,50.,od,skala_2);

  //pelne_energie:
  TH1D* hist_ph_pe = new TH1D("ph_pe","ph_pe;Liczba fotoelektronow;Liczba przypadkow",50.,0.,12000.); //1.- ile slupkow w histogramie, 2. i 3.- wymiary na OX
  TH1D* hist_ph_raw = new TH1D("ph_raw","ph_raw;Energia [ADC];Liczba przypadkow",50.,0.,800000.);
  TH1D* hist_ph_siglin = new TH1D("ph_siglin","ph_siglin;Energia [siglin];Liczba przypadkow",50.,0.,900000.);
  TH1D* hist_ph_sigcor = new TH1D("ph_sigcor","ph_sigcor;Energia [sigcor];Liczba przypadkow",50.,0.,900000.);

  //niskie_energie:
  TH1D* hist_ph_peLOW_CC = new TH1D("ph_pe","ph_pe;Liczba fotoelektronow;Liczba przypadkow CC",50.,0.,30.);//3000.); 
  TH1D* hist_ph_rawLOW_CC = new TH1D("ph_raw","ph_raw;Energia [ADC];Liczba przypadkow CC",50.,0.,3500.);//200000.);
  TH1D* hist_ph_siglinLOW_CC = new TH1D("ph_siglin","ph_siglin;Energia [siglin];Liczba przypadkow CC",50.,0.,3500.);//200000.);
  TH1D* hist_ph_sigcorLOW_CC = new TH1D("ph_sigcor","ph_sigcor;Energia [sigcor];Liczba przypadkow CC",50.,0.,3500.);//200000.);
  TH1D* hist_ph_peLOW_NC = new TH1D("ph_pe","ph_pe;Liczba fotoelektronow;Liczba przypadkow NC",50.,0.,30.);//3000.); 
  TH1D* hist_ph_rawLOW_NC = new TH1D("ph_raw","ph_raw;Energia [ADC];Liczba przypadkow NC",50.,0.,3500.);//200000.);
  TH1D* hist_ph_siglinLOW_NC = new TH1D("ph_siglin","ph_siglin;Energia [siglin];Liczba przypadkow NC",50.,0.,3500.);//200000.);
  TH1D* hist_ph_sigcorLOW_NC = new TH1D("ph_sigcor","ph_sigcor;Energia [sigcor];Liczba przypadkow NC",50.,0.,3500.);//200000.);

  //hist_NC_CC:
  TH1D* sigcor_NC = new TH1D ("sigcor NC","sigcor NC;Energia [sigcor];Liczba przypadkow",50.,0.,900000);
  TH1D* sigcor_CC = new TH1D ("sigcor CC","sigcor CC;Energia [sigcor];Liczba przypadkow",50.,0.,900000);
  TH1D* sigcor_UNKNOWN = new TH1D ("sigcor UNKNOWN","sigcor UNKNOWN;Energia [sigcor];Liczba przypadkow",50.,0.,900000);
  TH1D* sigcor_Suma = new TH1D ("sigcor Suma","sigcor Suma;Energia [sigcor];Liczba przypadkow",50.,0.,900000);

  //hist_ZERO_JEDEN:
  TH1D* sigcor_ZERO = new TH1D ("NC","NC;Energia [GeV];Liczba przypadkow",50.,0.,10.);//0
  TH1D* sigcor_JEDEN = new TH1D ("CC","CC;Energia [GeV];Liczba przypadkow",50.,0.,10.);//1

  //Pur_Eff:
  //zmienne niezbedne do histogramow Pur i Eff
  TH1F *NCall = new TH1F("NCall","NCall;Energia [GeV]; wszystkie NC",50,0,10);
  TH1F *CCall = new TH1F("CCall","CCall;Energia [GeV]; wszystkie CC",50,0,10);
  TH1F *CCjakoCC = new TH1F("CCjakoCC","CCjakoCC;Energia [GeV];CC jako CC",50,0,10); //pousuwac opisy osi
  TH1F *CCjakoNC = new TH1F("CCjakoNC","CCjakoNC;Energia [GeV];CC jako NC",50,0,10);
  TH1F *NCjakoNC = new TH1F("NCjakoNC","NCjakoNC;Energia [GeV];NC jako NC",50,0,10);
  TH1F *NCjakoCC = new TH1F("NCjakoCC","NCjakoCC;Energia [GeV];NC jako CC",50,0,10);
  CCjakoCC->Sumw2();
  CCjakoNC->Sumw2();
  NCjakoCC->Sumw2();
  NCjakoNC->Sumw2();
  NCall->Sumw2();
  CCall->Sumw2();

  //hist_STP:
  TH1D* stp_PH = new TH1D ("STP.PH.SIGCOR","STP.PH.SIGCOR;liczba plaszczyzn;Liczba przypadkow",20.,0.,20);
  //TH1D* stp_PH1 = new TH1D ("STP.PH1.SIGCOR","STP.PH1.SIGCOR;Energia [sigcor];Liczba przypadkow",50.,0.,650000);

  //galaz_MC:
  TH1D* mc_P4NEU_NC = new TH1D ("MC.P4NEU","MC.P4NEU;Energia [GeV];Liczba przypadkow",50.,0.,20);
  TH1D* mc_SUMA_NC = new TH1D ("SUMA","SUMA;Energia [GeV];Liczba przypadkow",50.,0.,20);
  TH1D* mc_P4NEU_CC = new TH1D ("MC.P4NEU","MC.P4NEU;Energia [GeV];Liczba przypadkow",50.,0.,20);
  TH1D* mc_SUMA_CC = new TH1D ("SUMA","SUMA;Energia [GeV];Liczba przypadkow",50.,0.,20);

  //profile_histogram:
  TProfile* profil = new TProfile ("evthdr.ph.sigcor VS mc.p4shw","evthdr.ph.sigcor VS mc.p4shw;evthdr.ph.sigcor;mc.p4shw",50.,0.,900000); //20 na koncu

  //hist_dlugosc:
  TH1D* dlugosc = new TH1D ("Dlugosc","Dlugosc;Dlugosc;Liczba przypadkow",50.,0.,900000);

  //////////////////////////// koniec deklaracji zmiennych

  cout<<endl<<"##############################"<<endl<<"Calkowita liczba przypadkow="<<nentries<<endl<<"Obliczenia... prosze czekac..."<<endl;










  //OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO
  ////////////////////////////////////////////////
  ////////////////////////////////////////////////
  /////////////////////START - GLOWNA PETLA PROGRAMU!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  // Loop over all tree entries
  ////////////////////////////////////////////////
  ////////////////////////////////////////////////
  //petla po wszystkich przypadkach:
  ////////////////////////////////////////////////
  //OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO
 







  for (Int_t jentry = 0; jentry < nentries; jentry++) 
  {
    Long64_t ientry = t.LoadTree(jentry);
    if (ientry < 0) break;
    nb = t.fChain->GetEntry(jentry);   nbytes += nb;
    int iaction = t.mc_iaction[0]; //= 1 dla CC; =0 dla NC - oznaczenie typu przypadku (tylko dla MC!)
    float enGeV = (double)t.evthdr_ph_sigcor*((double)0.0000834864); //energia podana w GeVach (ten przelicznik wziety z profile_histogram)
    int czastka=0, czastka1=0; //=1 dla CC; =0 dla NC; =-1 dla unknown; wspaznik jakiego typu dany przypadek PO selekcji 
    //if (jentry==10 || jentry==11 || jentry==12 || jentry==13) cout<<endl<<"planeall.beg: "<<t.evt_stp[3]<<endl;

    ///////histogramy wstepne://///////////////////
    
    //$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
    //hist_ZERO_JEDEN:
    if (t.evthdr_ph_sigcor>0 && GeV==1)
    {
      if (iaction==1) {jeden++; sigcor_JEDEN -> Fill(enGeV);} //CC
      if (iaction==0) {zero++; sigcor_ZERO -> Fill(enGeV);} //NC
    }
    //$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
/*    //pelne_energie:
    //if (iaction==0 )      //1 dla CC, 0 dla NC
    if (t.evthdr_ph_pe > 0)
    hist_ph_pe -> Fill(t.evthdr_ph_pe);
    if (t.evthdr_ph_raw > 0)
    hist_ph_raw -> Fill(t.evthdr_ph_raw);
    if (t.evthdr_ph_siglin > 0)
    hist_ph_siglin -> Fill(t.evthdr_ph_siglin);
    if (t.evthdr_ph_sigcor > 0)
    hist_ph_sigcor -> Fill(t.evthdr_ph_sigcor);
    //$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
    //niskie_energie:
    if (iaction==0) //NC
    {
      if (t.evthdr_ph_pe > 0)
      hist_ph_peLOW_NC -> Fill(t.evthdr_ph_pe);
      if (t.evthdr_ph_raw > 0)
      hist_ph_rawLOW_NC -> Fill(t.evthdr_ph_raw);
      if (t.evthdr_ph_siglin > 0)
      hist_ph_siglinLOW_NC -> Fill(t.evthdr_ph_siglin);
      if (t.evthdr_ph_sigcor > 0)
      hist_ph_sigcorLOW_NC -> Fill(t.evthdr_ph_sigcor);
    }
    if (iaction==1) //CC
    {
      if (t.evthdr_ph_pe > 0)
      hist_ph_peLOW_CC -> Fill(t.evthdr_ph_pe);
      if (t.evthdr_ph_raw > 0)
      hist_ph_rawLOW_CC -> Fill(t.evthdr_ph_raw);
      if (t.evthdr_ph_siglin > 0)
      hist_ph_siglinLOW_CC -> Fill(t.evthdr_ph_siglin);
      if (t.evthdr_ph_sigcor > 0)
      hist_ph_sigcorLOW_CC -> Fill(t.evthdr_ph_sigcor);
    }
    //$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
 */   //hist_STP:
    //zmienne pomocnicze do petli
/*    double suma1=0, suma2=0;
    for (int j=t.evthdr_planeall_beg; j<t.evthdr_planeall_end; j++)  //plane beg - plane.end
    { suma1=0; //suma2=0;
      //cout<<" "<<t.evthdr_nstrip;
      //if (jentry==10)
      for (int i=0;i<t.evthdr_nstrip;i++)
      {
         if (t.stp_plane[i]==j) suma1=suma1+t.stp_ph0_pe[i]+t.stp_ph1_pe[i];
         //cout<<" "<<t.stp_plane[i];
         //stp_PH -> Fill(suma1);
         //if (t.stp_ph1_sigcor[i]>0) suma2+=t.stp_ph1_sigcor[i];
         //stp_PH1 -> Fill(suma2);
      }
      if (suma1>3) suma2++;// cout<<"T";}
    } //*/
    //stp_PH -> Fill(suma2);  //pierwotnie suma1
 /*   //$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
    //galaz_MC:
    double Shower=0, suma3=0;
    for (int i=0; i<4; i++) Shower+=t.mc_p4shw[0][i];
    //cout<<"Shower= "<<Shower<<endl;
    //if (t.mc_Enu>0) mc_P4NEU->Fill(t.mc_Enu);
    if (t.mc_p4neu[0][3]>0 && iaction==0) mc_P4NEU_NC->Fill(t.mc_p4neu[0][3]);
    if (t.mc_p4neu[0][3]>0 && iaction==1) mc_P4NEU_CC->Fill(t.mc_p4neu[0][3]);
    suma3=t.mc_p4shw[0][3]+t.mc_p4mu1[0][3]+t.mc_p4mu2[0][3]+t.mc_p4el1[0][3]+t.mc_p4el2[0][3]+t.mc_p4tau[0][3];
    if (iaction==0) mc_SUMA_NC->Fill(suma3);
    if (iaction==1) mc_SUMA_CC->Fill(suma3);
    //+t.p4mu1[0][0]+t.p4mu2[0][0]+t.p4el1[0][0]+t.p4el2[0][0]+t.p4tau[0][0]);
    //if (t.mc_p4shw[0][0]>0) p4SHW->Fill(t.mc_p4shw[0][0]);
    //if (t.mc_tphv[0]>0) mc_TPHV->Fill(t.mc_tphv[0]);
    //if (t.mc_tphu[1]>0) mc_TPHU->Fill(t.mc_tphu[1]);
    //$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
    //profile_histogram:
    if (iaction==0) profil->Fill(t.evthdr_ph_sigcor,t.mc_p4shw[0][3]);
    //$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
    //koniec histogramow wstepnych
    */   

    //////////////////////////////////////////////////////////////////////////////////
    //rozpoczyna sie selekcja

    //obliczanie maksimow roznych wielkosci - na potrzeby zmiennych zlozonych
    double maks1=0, maks2=0, maks3=0, maks4=0;
    for (int i=0; i<=t.evthdr_ntrack; i++)
    {
      if (t.trk_range[i]>maks1) maks1=t.trk_range[i]; //zasieg wiodacego toru
      if (t.trk_ph_sigcor[i]>maks2) maks2=t.trk_ph_sigcor[i]; //w trackPHfraction - sygnal wiodacego toru. Calosc trkMaxPh
    }
    ////////////////
    for (int i=0; i<=t.evthdr_nshower; i++)
    {
      if (t.shw_ph_sigcor[i]>maks3) maks3=t.shw_ph_sigcor[i]; //shwMaxPh
      if (t.shw_ndigit[i]>maks4) maks4=t.shw_ndigit[i];
    } 
    if (maks4==0) maks4=1; //zabezbieczenie przed dzieleniem przez 0 - prawie nie gra roznicy
  
    //diffTrackShowerLengthsU - roznica miedzy dlugoscia wiodacego toru, a dlugoscia wiodacej kaskady U
    float maxShwRan=0;
    for(int i=0;i<t.evthdr_nshower;i++)
    {
      float ShwRanCur = t.shw_plane_endu[i]-t.shw_plane_begu[i];
      if(ShwRanCur>maxShwRan) maxShwRan = ShwRanCur;
    }
    
    float maxTrkRan=0;
    for(int i=0;i<t.evthdr_ntrack;i++)
    {
      float TrkRanCur = t.trk_plane_endu[i]-t.trk_plane_begu[i];
      if(TrkRanCur>maxTrkRan) maxTrkRan = TrkRanCur;
    }
  
    //diffTrackShowerLengthsV - roznica miedzy dlugoscia wiodacego toru, a dlugoscia wiodacej kaskady V
    float maxShwRanV=0;
    for(int i=0;i<t.evthdr_nshower;i++)
    {
      float ShwRanCurV = t.shw_plane_endv[i]-t.shw_plane_begv[i];
      if(ShwRanCurV>maxShwRanV) maxShwRanV = ShwRanCurV;
    }
   
    float maxTrkRanV=0;
    for(int i=0;i<t.evthdr_ntrack;i++)
    {
      float TrkRanCurV = t.trk_plane_endv[i]-t.trk_plane_begv[i];
      if(TrkRanCurV>maxTrkRanV) maxTrkRanV = TrkRanCurV;
    }
  
    //dlugosc przypadku
    double dlug=t.evthdr_plane_end-t.evthdr_plane_beg;
 
    //stosunek warosci sygnalu w roznych miejscach toru (liczba procentow zmienna)
    double pierwszy=0, drugi=0, pierwszy_prim=0;
    for (int i=0;i<t.evthdr_nstrip;i++)
    {
       if ( (t.stp_plane[i]>t.evthdr_plane_beg) && (t.stp_plane[i]<=t.evthdr_plane_beg+dlug*10./100.) )
       {
          pierwszy+=t.stp_ph0_sigcor[i]+t.stp_ph1_sigcor[i];
       }
       if ( (t.stp_plane[i]>t.evthdr_plane_beg) && (t.stp_plane[i]<=t.evthdr_plane_beg+dlug*25./100.) )
       {
          pierwszy_prim+=t.stp_ph0_sigcor[i]+t.stp_ph1_sigcor[i];
       }
       if ( (t.stp_plane[i]<t.evthdr_plane_end) && (t.stp_plane[i]>=t.evthdr_plane_end-dlug*33./100.) )
       {
          drugi+=t.stp_ph0_sigcor[i]+t.stp_ph1_sigcor[i];
       }
    }  


    if (t.evthdr_ph_sigcor>0 && TEST2==1) //histogram testowy TEST2 - rozklad dowolnej zmiennej przed selekcja
    {
       double zmienna_2=enGeV; 
       switch (zmiennaTEST)  //jakie zmienne bierzemy pod uwage przy N-1
       {
         case 28: zmienna_2=pierwszy_prim/maks3; break;
         case 3: zmienna_2=t.evthdr_plane_n; break;
         case 18: zmienna_2=maks4; break;
         case 26: zmienna_2=pierwszy/maks3; break;
         case 7: zmienna_2=t.evthdr_ph_sigcor/t.evthdr_planeall_n; break;
         case 2: zmienna_2=t.evthdr_ph_pe; break;
         case 13: zmienna_2=maxTrkRanV-maxShwRanV; break;
         case 20: zmienna_2=maks3; break;
         case 8: zmienna_2=maks1/t.evthdr_ph_sigcor; break;
         case 6: zmienna_2=t.evthdr_nstrip; break;
         case 0: zmienna_2=enGeV; break;
         case 91: zmienna_2=t.evthdr_ph_sigcor; break;
         case 92: zmienna_2=t.evthdr_ph_pe; break;
         case 93: zmienna_2=t.evthdr_nstrip; break;
         case 94: zmienna_2=t.evthdr_planeall_n; break;
         case 95: zmienna_2=t.evthdr_plane_n; break;
         case 96: zmienna_2=t.evthdr_ntrack; break;
         case 97: zmienna_2=t.evthdr_nshower; break;
       }

       if (iaction==0) TEST_NC2 -> Fill(zmienna_2);
       if (iaction==1) TEST_CC2 -> Fill(zmienna_2);
    }


    ////////////////////////////////////////////////
    //SELEKCJA WŁAŚCIWA:
    ////////////////////////////////////////////////
    czastka=-1; //czyli jaka, wg algorytmu selekcji, to jest czastka
    if (t.evthdr_ph_sigcor>0 && (PurEff==1 || N_1==1 || TEST1==1)) //warunek na niezerowa energie
    {
      if (iaction==1) ileCCwszystkich++; //oblicza ile jest przypadkow CC i NC dla sigcor>0
      if (iaction==0) ileNCwszystkich++;

      double cieciaXC=enGeV;

      switch (zmiennaN_1)  //jakie zmienne bierzemy pod uwage przy N-1
      {
        case 28: cieciaXC=pierwszy_prim/maks3; break;
        case 3: cieciaXC=t.evthdr_plane_n; break;
        case 18: cieciaXC=maks4; break;
        case 26: cieciaXC=pierwszy/maks3; break;
        case 7: cieciaXC=t.evthdr_ph_sigcor/t.evthdr_planeall_n; break;
        case 2: cieciaXC=t.evthdr_ph_pe; break;
        case 13: cieciaXC=maxTrkRanV-maxShwRanV; break;
        case 20: cieciaXC=maks3; break;
        case 8: cieciaXC=maks1/t.evthdr_ph_sigcor; break;
        case 6: cieciaXC=t.evthdr_nstrip; break;
      }

      //ciecia:
      //uwaga! zmienne licznik_wewCC oraz NC podaja liczbe KANDYDATOW na CC lub NC
      if (zmienna28==1 && pierwszy_prim/maks3>=0.52) //28
      {
         sigcor_CC -> Fill(t.evthdr_ph_sigcor); licznik_wewCC++; czastka=1;
         sigcor_Suma -> Fill(t.evthdr_ph_sigcor); if (iaction==CC) {N1_cieciaCC->Fill(cieciaXC); ileN_CC++;}
         ileCCpo28++; goto qqq; //continue; //break;
      }
      else
      if (zmienna28==1 && pierwszy_prim/maks3<=0.14) //zmienna 28 (pierwszy dla 25%)
      {
         sigcor_NC -> Fill(t.evthdr_ph_sigcor); licznik_wewNC++; czastka=0;
         sigcor_Suma -> Fill(t.evthdr_ph_sigcor); if (iaction==NC) {N1_cieciaNC->Fill(cieciaXC); ileN_NC++;}
         ileNCpo28++; goto qqq;// continue; // break;
      }
      else
  
      if (zmienna3==1 && t.evthdr_plane_n<=10) //3
      {
         sigcor_NC -> Fill(t.evthdr_ph_sigcor); licznik_wewNC++; czastka=0;
         sigcor_Suma -> Fill(t.evthdr_ph_sigcor); if (iaction==NC) {N1_cieciaNC->Fill(cieciaXC); ileN_NC++;}
         ileNCpo3++; goto qqq; //continue; // break;
      }
      else
      if (zmienna3==1 && t.evthdr_plane_n>=36) //zmienna 3
      {
         sigcor_CC -> Fill(t.evthdr_ph_sigcor); licznik_wewCC++; czastka=1;
         sigcor_Suma -> Fill(t.evthdr_ph_sigcor); if (iaction==CC) {N1_cieciaCC->Fill(cieciaXC); ileN_CC++;}
         ileCCpo3++; goto qqq; //continue; // break;
      }
      else
  
/*      if (zmienna18==1 && maks4<=60) //18
      {
         sigcor_CC -> Fill(t.evthdr_ph_sigcor); licznik_wewCC++; czastka=1;
         sigcor_Suma -> Fill(t.evthdr_ph_sigcor); if (iaction==CC) {N1_cieciaCC->Fill(cieciaXC); ileN_CC++;}
         ileCCpo18++; goto qqq; //continue; // break;
      }
      else
      if (zmienna18==1 && maks4>=400) //zmienna 18
      {
         sigcor_NC -> Fill(t.evthdr_ph_sigcor); licznik_wewNC++; czastka=0;
         sigcor_Suma -> Fill(t.evthdr_ph_sigcor); if (iaction==NC) {N1_cieciaNC->Fill(cieciaXC); ileN_NC++;}
         ileNCpo18++; goto qqq; //continue; // break;
      }
      else  //
  
      if (zmienna26==1 && pierwszy/maks3>=0.29) //zmienna 26 - uwaga na zmiany "pierwszy"!!!
      {
         sigcor_CC -> Fill(t.evthdr_ph_sigcor); licznik_wewCC++; czastka=1;
         sigcor_Suma -> Fill(t.evthdr_ph_sigcor); if (iaction==CC) {N1_cieciaCC->Fill(cieciaXC); ileN_CC++;}
         ileCCpo26++; goto qqq; //continue; // break;
      }
      else  
  */
      if (zmienna7==1 && t.evthdr_ph_sigcor/t.evthdr_planeall_n<=3400) //zmienna 7
      {
         sigcor_CC -> Fill(t.evthdr_ph_sigcor); licznik_wewCC++; czastka=1;
         sigcor_Suma -> Fill(t.evthdr_ph_sigcor); if (iaction==CC) {N1_cieciaCC->Fill(cieciaXC); ileN_CC++;}
         ileCCpo7bis++; goto qqq; //continue; // break;
      }
      else
      if (zmienna7==1 && t.evthdr_ph_sigcor/t.evthdr_planeall_n>=7800) //zmienna 7
      {
         sigcor_CC -> Fill(t.evthdr_ph_sigcor); licznik_wewCC++; czastka=1;
         sigcor_Suma -> Fill(t.evthdr_ph_sigcor); if (iaction==CC) {N1_cieciaCC->Fill(cieciaXC); ileN_CC++;}
         ileCCpo7bis++; goto qqq; //continue; // break;
      }
      else
/*  
      if (zmienna2==1 && t.evthdr_ph_pe<=500) //zmienna 2
      {
         sigcor_NC -> Fill(t.evthdr_ph_sigcor); licznik_wewNC++; czastka=0;
         sigcor_Suma -> Fill(t.evthdr_ph_sigcor); if (iaction==NC) {N1_cieciaNC->Fill(cieciaXC); ileN_NC++;}
         ileNCpo2++; goto qqq; //continue; // break;
      } 
      else  //
  
      if (zmienna13==1 && maxTrkRanV-maxShwRanV>=9) //zmienna 13
      {
         sigcor_CC -> Fill(t.evthdr_ph_sigcor); licznik_wewCC++; czastka=1;
         sigcor_Suma -> Fill(t.evthdr_ph_sigcor); if (iaction==CC) {N1_cieciaCC->Fill(cieciaXC); ileN_CC++;}
         ileCCpo13++; goto qqq; //continue; // break;
      }
      else
  */
      if (zmienna20==1 && (maks3>=94000 || maks3<=50000)) //20
      {
         sigcor_CC -> Fill(t.evthdr_ph_sigcor); licznik_wewCC++; czastka=1;
         sigcor_Suma -> Fill(t.evthdr_ph_sigcor); if (iaction==CC) {N1_cieciaCC->Fill(cieciaXC); ileN_CC++;}
         ileCCpo20++; goto qqq; //continue; // break;
      }
      else
      if (zmienna20==1 && maks3>=60000 && maks3<=70000) //zmienna 20
      {
         sigcor_NC -> Fill(t.evthdr_ph_sigcor); licznik_wewNC++; czastka=0;
         sigcor_Suma -> Fill(t.evthdr_ph_sigcor); if (iaction==NC) {N1_cieciaNC->Fill(cieciaXC); ileN_NC++;}
         ileNCpo20++; goto qqq; //continue; // break;
      }
      else 
      if (zmienna20==1 && maks3>=86000 && maks3<94000) //zmienna 20
      {
         sigcor_NC -> Fill(t.evthdr_ph_sigcor); licznik_wewNC++; czastka=0;
         sigcor_Suma -> Fill(t.evthdr_ph_sigcor); if (iaction==NC) {N1_cieciaNC->Fill(cieciaXC); ileN_NC++;}
         ileNCpo20++; goto qqq; //continue; // break;
      }
      else 
      if (zmienna20==1 && maks3>70000 && maks3<86000) //20
      {
         sigcor_CC -> Fill(t.evthdr_ph_sigcor); licznik_wewCC++; czastka=1;
         sigcor_Suma -> Fill(t.evthdr_ph_sigcor); if (iaction==CC) {N1_cieciaCC->Fill(cieciaXC); ileN_CC++;}
         ileCCpo20++; goto qqq; //continue; // break;
      }
      else
/*  
      if (zmienna8==1 && maks1/t.evthdr_ph_sigcor>=0.015) //8
      {
         sigcor_CC -> Fill(t.evthdr_ph_sigcor); licznik_wewCC++; czastka=1;
         sigcor_Suma -> Fill(t.evthdr_ph_sigcor); if (iaction==CC) {N1_cieciaCC->Fill(cieciaXC); ileN_CC++;}
         ileCCpo8++; goto qqq; //continue; // break;
      }
      else
      if (zmienna8==1 && maks1/t.evthdr_ph_sigcor<=0.0045) //zmienna 8
      {
         sigcor_NC -> Fill(t.evthdr_ph_sigcor); licznik_wewNC++; czastka=0;
         sigcor_Suma -> Fill(t.evthdr_ph_sigcor); if (iaction==NC) {N1_cieciaNC->Fill(cieciaXC); ileN_NC++;}
         ileNCpo8++; goto qqq; //continue; // break;
      }
      else
*/  
      if (zmienna6==1 && t.evthdr_nstrip>=62) //6
      {
         sigcor_NC -> Fill(t.evthdr_ph_sigcor); licznik_wewNC++; czastka=0;
         sigcor_Suma -> Fill(t.evthdr_ph_sigcor); if (iaction==NC) {N1_cieciaNC->Fill(cieciaXC); ileN_NC++;}
         ileNCpo20++; goto qqq; //continue; // break;
      }   // */
/*      else 
      //to, co pozostalo po selekcji:
      {
         sigcor_CC -> Fill(t.evthdr_ph_sigcor); licznik_wewCC++; czastka=1;
         sigcor_Suma -> Fill(t.evthdr_ph_sigcor);
         ileCCelse++; goto qqq; //continue; // break;
      }  // */
      else {czastka=-1; goto qqq;} //unknown

    } //<- koniec ciec


    //tu jestem juz po cieciach na zmiennych. Mam KANDYDATOW na CC/NC
    ///////////////////////////////////////////////////////////////////////////////

    qqq: //po kazdym cieciu skacze tutaj

    if (TEST1==1)
    {
    //histogram testowy nr 1 - po cieciach:
    double zmienna=enGeV;
    switch (zmiennaTEST)  //jakie zmienne bierzemy pod uwage przy N-1
    {
      case 28: zmienna=pierwszy_prim/maks3; break;
      case 3: zmienna=t.evthdr_plane_n; break;
      case 18: zmienna=maks4; break;
      case 26: zmienna=pierwszy/maks3; break;
      case 7: zmienna=t.evthdr_ph_sigcor/t.evthdr_planeall_n; break;
      case 2: zmienna=t.evthdr_ph_pe; break;
      case 13: zmienna=maxTrkRanV-maxShwRanV; break;
      case 20: zmienna=maks3; break;
      case 8: zmienna=maks1/t.evthdr_ph_sigcor; break;
      case 6: zmienna=t.evthdr_nstrip; break;
      case 0: zmienna_2=enGeV; break;
      case 91: zmienna_2=t.evthdr_ph_sigcor; break;
      case 92: zmienna_2=t.evthdr_ph_pe; break;
      case 93: zmienna_2=t.evthdr_nstrip; break;
      case 94: zmienna_2=t.evthdr_planeall_n; break;
      case 95: zmienna_2=t.evthdr_plane_n; break;
      case 96: zmienna_2=t.evthdr_ntrack; break;
      case 97: zmienna_2=t.evthdr_nshower; break;
    }

    if (czastka==0) {TEST_NC -> Fill(zmienna);}
    if (czastka==1) {TEST_CC -> Fill(zmienna);}
    }

    ///////////////////////////////////////////Zapasowa (stara) selekcja (przydatne do generowania histogramow dodatkowych
/*
    //double dlug=t.evthdr_plane_end-t.evthdr_plane_beg;

    //   if (pierwszy_prim/maks3>=0.63 || pierwszy_prim/maks3<=0.13 || t.evthdr_plane_n<=9 || t.evthdr_plane_n>=37 || maks4<=60 || maks4>=250 || pierwszy/maks3>=0.28 || t.evthdr_ph_sigcor/t.evthdr_planeall_n<=2400 || t.evthdr_ph_pe<=500 || maxTrkRanV-maxShwRanV>=5 || maks3>=85000 || maks3<=48000 || maks1/t.evthdr_ph_sigcor>=0.015 || maks1/t.evthdr_ph_sigcor<=0.005 || (t.evthdr_nstrip>=95 && t.evthdr_nstrip<=102))   
       
       czastka=-1;
        
       //wstawic dodatkowe warunki w IFach!!!
       if (pierwszy_prim/maks3>=0.63) {czastka=1; goto q1;} else //zmienna 28
       if (pierwszy_prim/maks3<=0.13) {czastka=0; goto q1;} else //zmienna 28 (pierwszy dla 25%)
       if (t.evthdr_plane_n<=9) {czastka=0; goto q1;} else //zmienna 3
       if (t.evthdr_plane_n>=37) {czastka=1; goto q1;} else //zmienna 3
         if (maks4<=60) {czastka=1; goto q1;} else //zmienna 18
         if (maks4>=250) {czastka=0; goto q1;} else //zmienna 18
       if (pierwszy/maks3>=0.28) {czastka=1; goto q1;} else //zmienna 26 - uwaga na zmiany "pierwszy"!!!
       if (t.evthdr_ph_sigcor/t.evthdr_planeall_n<=2400) {czastka=1; goto q1;} else //zmienna 7 - ciecie do przemyslenia!
         if (t.evthdr_ph_pe<=500) {czastka=0; goto q1;} else //zmienna 2
       if (maxTrkRanV-maxShwRanV>=5) {czastka=1; goto q1;} else //zmienna 13
       if (maks3>=85000) {czastka=1; goto q1;} else //zmienna 20
       if (maks3<=48000) {czastka=0; goto q1;} else //zmienna 20
       if (maks1/t.evthdr_ph_sigcor>=0.015) {czastka=1; goto q1;} else //zmienna 8
       if (maks1/t.evthdr_ph_sigcor<=0.005) {czastka=0; goto q1;} else //zmienna 8 
       if (t.evthdr_nstrip>=95 && t.evthdr_nstrip<=102) {czastka=0; goto q1;} else {czastka=-1;} //zmienna 6  

       q1:
       //if (t.evthdr_ph_sigcor>0)
       {
          //zmienne ileN_NC i CC podaja ile jest rzeczywistych przypadkow NC i CC posrod PELNEJ puli wyselekcjonowanych
          //double zmienna=t.evthdr_plane_n;///do histogramu TEST nr 1 
          double cieciaXC=enGeV;
          switch (zmiennaN_1)  //jakie zmienne bierzemy pod uwage przy N-1
          {
            case 28: cieciaXC=pierwszy_prim/maks3; break;
            case 3: cieciaXC=t.evthdr_plane_n; break;
            case 18: cieciaXC=maks4; break;
            case 26: cieciaXC=pierwszy/maks3; break;
            case 7: cieciaXC=t.evthdr_ph_sigcor/t.evthdr_planeall_n; break;
            case 2: cieciaXC=t.evthdr_ph_pe; break;
            case 13: cieciaXC=maxTrkRanV-maxShwRanV; break;
            case 20: cieciaXC=maks3; break;
            case 8: cieciaXC=maks1/t.evthdr_ph_sigcor; break;
            case 6: cieciaXC=t.evthdr_nstrip; break;
          }

          if (czastka==0 && iaction==NC) {/*TEST_NC -> Fill(zmienna); N1_cieciaNC->Fill(cieciaXC); ileN_NC++;}
          if (czastka==1 && iaction==CC) {/*TEST_CC -> Fill(zmienna); N1_cieciaCC->Fill(cieciaXC); ileN_CC++;}
          //cout<<"Snarls: "<<t.detsim_snarlDigits<<" Nr przypadku: "<<jentry<<endl;
       }       
    }

    //q2:     */

    //////////////////////////////////////////////
    //czystosc i efektrywnosc - liczenie wartosci
    //////////////////////////////////////////////

    if (/*t.evthdr_ph_sigcor>0 &&*/ PurEff==1 ) //tyczy sie tylko Eff (mianownik!!!)
    {
       if ( iaction == 1 )
       {
          iIleCC++; //jeden
          CCall->Fill(enGeV);
       }
       else 
       if ( iaction == 0 )
       {
          iIleNC++; //zero
          NCall->Fill(enGeV);
       }
    }    

    if (t.evthdr_ph_sigcor>0 && PurEff==1)
    {
       if ( iaction == 1 || iaction == 0 )
       {
          if (iaction == 1)
          { //CC
             if (czastka == 1)
             {
                itCCjakoCC++;
                CCjakoCC->Fill(enGeV);
             }
             else 
             if (czastka == 0) //==0 lub -1
             { //NC
                itCCjakoNC++;
                CCjakoNC->Fill(enGeV);
             }
          }  
          else  //czyli iaction==0
          {//tutaj NC
             if (czastka == 0)
             {
                itNCjakoNC++;
                NCjakoNC->Fill(enGeV);
             }
             else 
             if (czastka == 1) //==1 lub -1
             {
                itNCjakoCC++;
                NCjakoCC->Fill(enGeV);
             }
          }  
       }  
    } //petla dla Pur i Eff
    ///////////////////////////////////////////////

    
    ///////////////////////////////histogram N-1

    if (t.evthdr_ph_sigcor>0 && N_1==1)
    {

       //histogram dla wszystkich przypadkow (niezbedne do wykresu N-1)

       double wszystkoXC=enGeV;
       switch (zmiennaN_1)  //jakie zmienne bierzemy pod uwage przy N-1
       {
         case 28: wszystkoXC=pierwszy_prim/maks3; break;
         case 3: wszystkoXC=t.evthdr_plane_n; break;
         case 18: wszystkoXC=maks4; break;
         case 26: wszystkoXC=pierwszy/maks3; break;
         case 7: wszystkoXC=t.evthdr_ph_sigcor/t.evthdr_planeall_n; break;
         case 2: wszystkoXC=t.evthdr_ph_pe; break;
         case 13: wszystkoXC=maxTrkRanV-maxShwRanV; break;
         case 20: wszystkoXC=maks3; break;
         case 8: wszystkoXC=maks1/t.evthdr_ph_sigcor; break;
         case 6: wszystkoXC=t.evthdr_nstrip; break;
       }
       if (iaction==0) {N1_wszystkoNC->Fill(wszystkoXC); ileALL_NC++;}
       if (iaction==1) {N1_wszystkoCC->Fill(wszystkoXC); ileALL_CC++;}
    

//       if (/*pierwszy_prim/maks3>=0.63 || pierwszy_prim/maks3<=0.13 || t.evthdr_plane_n<=9 || t.evthdr_plane_n>=37 || maks4<=60 || maks4>=250 || pierwszy/maks3>=0.28 || t.evthdr_ph_sigcor/t.evthdr_planeall_n<=2400 || t.evthdr_ph_pe<=500 || maxTrkRanV-maxShwRanV>=5 || maks3>=85000 || maks3<=48000 || maks1/t.evthdr_ph_sigcor>=0.015 || maks1/t.evthdr_ph_sigcor<=0.005 || (t.evthdr_nstrip>=95 && t.evthdr_nstrip<=102))   


       czastka1=0; //0=NC; 1=CC

       if (zmienna28==1 && zmiennaN_1!=28 && pierwszy_prim/maks3>=0.52) {czastka1=1; goto q3;} else //zmienna 28
       if (zmienna28==1 && zmiennaN_1!=28 && pierwszy_prim/maks3<=0.14) {czastka1=0; goto q3;} else //zmienna 28 (pierwszy dla 25%)
       if (zmienna3==1 && zmiennaN_1!=3 && t.evthdr_plane_n<=10) {czastka1=0; goto q3;} else //zmienna 3
       if (zmienna3==1 && zmiennaN_1!=3 && t.evthdr_plane_n>=36) {czastka1=1; goto q3;} else //zmienna 3
//         if (zmienna18==1 && zmiennaN_1!=18 && maks4<=60) {czastka1=1; goto q3;} else //zmienna 18
//         if (zmienna18==1 && zmiennaN_1!=18 && maks4>=400) {czastka1=0; goto q3;} else //zmienna 18
//       if (zmienna26==1 && zmiennaN_1!=26 && pierwszy/maks3>=0.29) {czastka1=1; goto q3;} else //zmienna 26 - uwaga na zmiany "pierwszy"!!!
       if (zmienna7==1 && zmiennaN_1!=7 && t.evthdr_ph_sigcor/t.evthdr_planeall_n<=3400) {czastka1=1; goto q3;} else //zmienna 7
       if (zmienna7==1 && zmiennaN_1!=7 && t.evthdr_ph_sigcor/t.evthdr_planeall_n>=7800) {czastka1=1; goto q3;} else //zmienna 7
//         if (zmienna2==1 && zmiennaN_1!=2 && t.evthdr_ph_pe<=500) {czastka1=0; goto q3;} else //zmienna 2
//       if (zmienna13==1 && zmiennaN_1!=13 && maxTrkRanV-maxShwRanV>=9) {czastka1=1; goto q3;} else //zmienna 13
       if (zmienna20==1 && zmiennaN_1!=20 && (maks3>=94000 || maks3<=50000)) {czastka1=1; goto q3;} else //zmienna 20
       if (zmienna20==1 && zmiennaN_1!=20 && maks3>=60000 && maks3<=70000) {czastka1=0; goto q3;} else //zmienna 20
       if (zmienna20==1 && zmiennaN_1!=20 && maks3>=86000 && maks3<94000) {czastka1=0; goto q3;} else //zmienna 20
       if (zmienna20==1 && zmiennaN_1!=20 && maks3>70000 && maks3<86000) {czastka1=1; goto q3;} else //zmienna 20
//       if (zmienna8==1 && zmiennaN_1!=8 && maks1/t.evthdr_ph_sigcor>=0.015) {czastka1=1; goto q3;} else //zmienna 8
//       if (zmienna8==1 && zmiennaN_1!=8 && maks1/t.evthdr_ph_sigcor<=0.0045) {czastka1=0; goto q3;} else //zmienna 8 
       if (zmienna6==1 && zmiennaN_1!=6 && t.evthdr_nstrip>=62) {czastka1=0; goto q3;} //zmienna 6  
       else {czastka1=-1;}
 
       q3:
       {
         double ciecieXC=enGeV;
         switch (zmiennaN_1)  //jakie zmienne bierzemy pod uwage przy N-1
         {
           case 28: ciecieXC=pierwszy_prim/maks3; break;
           case 3: ciecieXC=t.evthdr_plane_n; break;
           case 18: ciecieXC=maks4; break;
           case 26: ciecieXC=pierwszy/maks3; break;
           case 7: ciecieXC=t.evthdr_ph_sigcor/t.evthdr_planeall_n; break;
           case 2: ciecieXC=t.evthdr_ph_pe; break;
           case 13: ciecieXC=maxTrkRanV-maxShwRanV; break;
           case 20: ciecieXC=maks3; break;
           case 8: ciecieXC=maks1/t.evthdr_ph_sigcor; break;
           case 6: ciecieXC=t.evthdr_nstrip; break;
         }

         if (czastka1==0 && iaction==NC) {N1_1ciecieNC->Fill(ciecieXC); ileN_1_NC++;}
         if (czastka1==1 && iaction==CC) {N1_1ciecieCC->Fill(ciecieXC); ileN_1_CC++;}
       }
    }  

 //   q4: 
   // if (t.evthdr_ph_sigcor<0) {} //aby oszukac "goto"

    ////////koniec N-1

    ////////// koniec TESTu

    //////////////////////////////////////////////////
    //"PO" - histogramy po poprzednich cieciach

    //selekcja ogolna
    if (t.evthdr_ph_sigcor>0 && H_PO==1)
    {
         double jakaPO=enGeV;
         switch (zmiennaPO)  //jakie zmienne bierzemy pod uwage przy N-1
         {
           case 28: jakaPO=pierwszy_prim/maks3; goto po1;
           case 3: jakaPO=t.evthdr_plane_n; 
                   if (zmienna28==1 && pierwszy_prim/maks3<0.52 && pierwszy_prim/maks3>0.14) //zmienna 28
                   goto po1;
           case 7: jakaPO=t.evthdr_ph_sigcor/t.evthdr_planeall_n; 
                   if (zmienna28==1 && pierwszy_prim/maks3<0.52 && pierwszy_prim/maks3>0.14) //zmienna 28
                   if (zmienna3==1 && t.evthdr_plane_n>10 && t.evthdr_plane_n<36) //zmienna 3
              //     if (zmienna18==1 && maks4>60 && maks4<400) //zmienna 18
              //     if (zmienna26==1 && pierwszy/maks3<0.29) //zmienna 26
                   goto po1;
           case 20: jakaPO=maks3; 
                    if (zmienna28==1 && pierwszy_prim/maks3<0.52 && pierwszy_prim/maks3>0.14) //zmienna 28
                    if (zmienna3==1 && t.evthdr_plane_n>10 && t.evthdr_plane_n<36) //zmienna 3
              //      if (zmienna18==1 && maks4>60 && maks4<400) //zmienna 18
              //      if (zmienna26==1 && pierwszy/maks3<0.29) //zmienna 26
                    if (zmienna7==1 && t.evthdr_ph_sigcor/t.evthdr_planeall_n>3400 && t.evthdr_ph_sigcor/t.evthdr_planeall_n<7800) //zmienna 7
              //      if (zmienna2==1 && t.evthdr_ph_pe>500) //zmienna 2
              //      if (zmienna13==1 && maxTrkRanV-maxShwRanV<9) //zmienna 13
                    goto po1;
           case 6: jakaPO=t.evthdr_nstrip; 
                   if (zmienna28==1 && pierwszy_prim/maks3<0.52 && pierwszy_prim/maks3>0.14) //zmienna 28
                   if (zmienna3==1 && t.evthdr_plane_n>10 && t.evthdr_plane_n<36) //zmienna 3
              //     if (zmienna18==1 && maks4>60 && maks4<400) //zmienna 18
              //     if (zmienna26==1 && pierwszy/maks3<0.29) //zmienna 26
                   if (zmienna7==1 && t.evthdr_ph_sigcor/t.evthdr_planeall_n>3400 && t.evthdr_ph_sigcor/t.evthdr_planeall_n<7800) //zmienna 7
              //     if (zmienna2==1 && t.evthdr_ph_pe>500) //zmienna 2
              //     if (zmienna13==1 && maxTrkRanV-maxShwRanV<9) //zmienna 13
                   if (zmienna20==1 && maks3<60000 && maks3>50000) //zmienna 20
              //     if (zmienna20==1 && (maks3<58000 || maks3>81000)) //zmienna 20
              //     if (zmienna8==1 && maks1/t.evthdr_ph_sigcor<0.015 && maks1/t.evthdr_ph_sigcor>0.0045) //zmienna 8 
                   goto po1;
           case 0: jakaPO=enGeV;
                   if (zmienna28==1 && pierwszy_prim/maks3<0.52 && pierwszy_prim/maks3>0.14) //zmienna 28
                   if (zmienna3==1 && t.evthdr_plane_n>10 && t.evthdr_plane_n<36) //zmienna 3
              //     if (zmienna18==1 && maks4>60 && maks4<400) //zmienna 18
              //     if (zmienna26==1 && pierwszy/maks3<0.29) //zmienna 26
                   if (zmienna7==1 && t.evthdr_ph_sigcor/t.evthdr_planeall_n>3400 && t.evthdr_ph_sigcor/t.evthdr_planeall_n<7800) //zmienna 7
              //     if (zmienna2==1 && t.evthdr_ph_pe>500) //zmienna 2
              //     if (zmienna13==1 && maxTrkRanV-maxShwRanV<9) //zmienna 13
                   if (zmienna20==1 && maks3<60000 && maks3>50000) //zmienna 20
              //     if (zmienna20==1 && (maks3<58000 || maks3>81000)) //zmienna 20
              //     if (zmienna8==1 && maks1/t.evthdr_ph_sigcor<0.015 && maks1/t.evthdr_ph_sigcor>0.0045) //zmienna 8 
                   if (zmienna6==1 && t.evthdr_nstrip<62) //zmienna 6 
                   goto po1;
           default: goto po2;
         }
         po1:
         {if (iaction==0) PO_NC->Fill(jakaPO);
         if (iaction==1) PO_CC->Fill(jakaPO);} 
         po2: 
         if (t.evthdr_ph_sigcor<0) {} 
    }

  } 

  //OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO
  ///////////////////////////////////////////////////////////////////////////
  ///////////////////////////////////////////////////////////////////////////
  ///////////////////////////////////////////////////////koniec glownej petli (po przypadkach)
  ///////////////////////////////////////////////////////////////////////////
  ///////////////////////////////////////////////////////////////////////////
  //OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO

  // czystosc i efektywnosc - wartosci
  if (PurEff==1)
  {
    float PurCC = (float) itCCjakoCC / (itCCjakoCC + itNCjakoCC);
    float EffCC = (float) itCCjakoCC / iIleCC;	
    float PurNC = (float) itNCjakoNC / (itNCjakoNC + itCCjakoNC);
    float EffNC = (float) itNCjakoNC / iIleNC;
  }

  //ileN_CC=ileN_CC+298;
  //ileN_1_CC=ileN_1_CC+298;
  /////////////////////////////////////////////////
  /////////////////////////////RAPORT:
  cout<<endl<<endl<<"###################################### Dane podstawowe:"<<endl;
  cout<<"ilosc CC mc (sigcor>0): "<<jeden<<" Ilosc NC mc: "<<zero<<endl;
  cout<<"ilosc CC po selek (licznik wewn) - czyli liczba KANDYDATOW: "<<licznik_wewCC<<" Ilosc NC po selek: "<<licznik_wewNC<<" Ilosc UNKNOWN: "<<licznik_wewUNKNOWN<<endl;
  cout<<"Liczba kandydatow: Ile N-1 CC: "<<ileN_1_CC<<" Ile N-1 NC: "<<ileN_1_NC<<" Ile N NC: "<<ileN_NC<<" Ile N CC: "<<ileN_CC<<" Ile wszystkich NC (mc): "<<ileALL_NC<<" Ile wszystkich CC (mc): "<<ileALL_CC<<endl;
  cout<<"Suma N-1: "<<ileN_1_CC+ileN_1_NC<<" Suma N: "<<ileN_NC+ileN_CC<<" Suma wszystkich: "<<ileALL_NC+ileALL_CC<<endl;
  cout<< ileN_1_CC+ileN_1_NC <<" (CC: "<<ileN_1_CC<<", NC: "<<ileN_1_NC<<")"<<endl;
  //cout<<"Licznik wewNC: "<<licznik_wewNC<<" Licznik wewCC: "<<licznik_wewCC<<endl;
  cout<<"###################################### Szczegolowe dane do Pur i Eff:"<<endl;
//  cout<<"CCjakoCC; CCjakoNC; NCjakoNC; NCjakoCC:"<<endl;
  cout << "CCjakoCC: "<< itCCjakoCC << " CCjakoNC: " << itCCjakoNC <<" NCjakoNC: ";
  cout << itNCjakoNC << " NCjakoCC: " << itNCjakoCC << endl;
  cout <<"IleCC; ileNC (mc) bez warunku na sigcor>0!: "<<iIleCC<<" ; "<<iIleNC<<endl;
  //cout << "Przypadkow CC mc: " << iIleCC << " NC mc: " << iIleNC << endl;
  if (PurEff==1)
  {
    cout << "################################### Pur i Eff:"<<endl;
    cout << "Czystosc CC " << PurCC << endl;
    cout << "Efektywnosc CC " << EffCC << endl;
    cout << "Czystosc NC " << PurNC << endl;
    cout << "Efektywnosc NC " << EffNC << endl<<endl<<endl;
  }
  cout << "################################### Dane po kazdym cieciu (oprocz 2 pierwszych, reszta to liczba KANDYDATOW):"<<endl;
  cout<<"ileNCwszystkich (mc i sigcor>0)="<<ileNCwszystkich<<" ileCCwszystkich (mc i sigcor>0)="<<ileCCwszystkich<<" ileNCpo28="<<ileNCpo28<<" ileCCpo28="<<ileCCpo28<<" ileNCpo3="<<ileNCpo3<<" ileCCpo3="<<ileCCpo3<<" ileNCpo18="<<ileNCpo18<<" ileCCpo18="<<ileCCpo18<<" ileCCpo26="<<ileCCpo26<<" ileCCpo7="<<ileCCpo7bis<<" ileNCpo2="<<ileNCpo2<<" ileCCpo13="<<ileCCpo13<<" ileCCpo20="<<ileCCpo20<<" ileNCpo20="<<ileNCpo20<<" ileCCpo8="<<ileCCpo8<<" ileNCpo8="<<ileNCpo8<<" ileNCpo6="<<ileNCpo6<<" ileCCelse="<<ileCCelse<<endl<<endl;  
  //cout<<"Dlug="<<dlug<<" Pierwszy="<<pierwszy<<" Drugi="<<drugi<<endl;

  ////////////////////////////////////////////////////////////
  //Rysowanie na ekranie:

  if (PurEff==1) Pur_Eff(CCjakoCC,CCjakoNC,NCjakoNC,NCjakoCC,CCall,NCall,itCCjakoCC,itCCjakoNC,itNCjakoNC,itNCjakoCC,iIleCC,iIleNC); //Histogramy Pur i Eff
  //pelne_energie(); //4 histogramy czterec zmiennych energetycznych
  //hist_NC_CC(); //selekcja NC/CC
  if (GeV==1) hist_ZERO_JEDEN(sigcor_JEDEN, sigcor_ZERO); //mc.iaction - rozklad energii NC/CC
  //hist_dlugosc(); //dlugosc przypadku
  //niskie_energie(); //histogramy dla niskich energii
  if (H_STP==1) hist_STP(stp_PH); //dane z galezi STP
  //galaz_MC(); //suma danych z galezi MC
  //profile_histogram(); //podaje zaleznosc miedzy sigcor a GeV
  if (TEST1==1) hist_TEST(); //test zmiennej NC/CC po cieciach
  if (TEST2==1) hist_TEST2(); //test zmiennej NC/CC bez ciec
  if (H_PO==1) HIST_PO(PO_NC,PO_CC); //histogram zmiennej po poprzednich cieciach
  if (N_1==1) Nmniej1(N1_wszystkoNC, N1_1ciecieNC, N1_cieciaNC, N1_wszystkoCC, N1_1ciecieCC, N1_cieciaCC); //N-1
}
///////////////////////////////////QNIEC