//    1. Bayesian fit (11) better than Pandel fit (7), by a factor specific
//for atmospheric neutrino hypothesis - (jkchi(11)-(jkchi(7)+29.67)>0). The
//equivalent for the old Bayes would be to just take those reconstructed
//upgoing.
//
//    2. ldirb(7)>100
//    3. abs(smootallphit(7))<0.3
//    4. abs(zenith(1)-zenith(7))<20  linefit - best pandel difference
//    5. ndirc(7)>7
//
//   Then I added -
//
//    nt/cuts $50 ndirc(7)>((100.*abs(smootallphit(7)))-10.)
//    nt/cuts $51 nch>((50.*abs(smootallphit(7)))+10.)
//
#include "../root/Common.h"

// mein Cuts sind natuerlich GLOBAL


TCut ZenZyl;
TCut LogENch;
TCut NLdirc;
TCut LCor;
TCut Lupdown;
TCut LTrCasc;
TCut MultiSmooth;
TCut PsiOne;
TCut PsiTwo;
TCut ZenZylHard;
TCut ZenAll;
TCut ZenMean;
TCut PsiMean;

void Init( char * option = "" ){
  cout << "Loading L3 File ...\n";
#ifndef plotter
  gROOT->LoadMacro("../root/plot.C");
#endif
  TFile * f = new TFile("/afs/ifh.de/user/s/sboeser/projects/ana2000"
                        "/data/R4_F3/root/L3.root",option);
  minset = (int)cexp; maxset = (int)coinc;
  f->cd(); f->ls();
  for ( int ds= minset; ds <= maxset; ds++){
    cout << "Getting " << (char*)Names[ds] << " from " 
         << f->Get((char*)Names[ds]) << "...\n";
    if ( Trees[ds] ) Trees[ds]->Delete();
    Trees[ds] = (TTree*)(f->Get((char*)Names[ds]));
    if ( ! Trees[ds] ){
      cout << "Trees " << Names[ds] << " not found\n ";
      cout.flush();
      break;
    }
  }
  cout << "Defining cuts !\n";
  // erstmal kommen MEINE cuts
  cout << " -> ZenZyl\n";
  ZenZyl.SetTitle("Zenith[10]>90");

  cout << " -> ZenZylHard\n";
  ZenZylHard.SetTitle("Zenith[10]>100");

  cout << " -> ZenAll\n";
  ZenAll.SetTitle("(Zenith[10]>90)&&(Zenith[6]>90)&&(Zenith[3]>90)&&(Zenith[1] > 90)");

  cout << " -> ZenMean\n";
  ZenMean.SetTitle("Mean.zen > 95");
	
  cout << " -> PsiMean\n";
  PsiMean.SetTitle("Mean.psi < .2");

  cout << " -> LogENch\n";
  LogENch.SetTitle("(Nch - 10*log10(Energy[7]))>0");

  cout << " -> NLdirc\n";
  NLdirc.SetTitle("Ndirc[10]*Ldirc[10] < 1000");

  cout << " -> LCor\n";
  LCor.SetTitle("(37.5*Ldirc[10]-Jkprob[5]) > 3000");

  cout << " -> Lupdown\n";
  Lupdown.SetTitle("Jkprob[10]-Jkchi[10] < 5");

  cout << " -> LTrCasc\n";
  LTrCasc.SetTitle("(Jkrchi[8]-Jkrchi[10]) >0 ");

  cout << " -> MultiSmooth\n";
  MultiSmooth.SetTitle("(Lvsu1/Nhitac1+abs(Szyldir1)+abs(Smootdirct[10]))<1.2");

  cout << " -> PsiOne\n";
  PsiOne.SetTitle(SolidAngle("Zenith[1]","Zenith[6]",
                    "Azimuth[1]","Azimuth[6]")+TString("<10"));
  cout << " -> PsiTwo\n";
  PsiTwo.SetTitle(SolidAngle("Zenith[10]","Zenith[3]",
                         "Azimuth[10]","Azimuth[3]")+TString("<20"));

  cout << " ... done\n";
}
 

void UnSelect (){
  for ( int ds = minset; ds <= maxset; ds++){
    cout << "Reseting Tree \n";
    Trees[ds]->SetEventList(0);
	}
}

void Select (TCut Sel){
  
  UnSelect();
	
	TH1F * hdummy = new TH1F("hdummy","Ich bin ein DUMMY",1,300,301);
  for ( int ds = minset; ds <= maxset; ds++){
    cout << "Selecting in " << (char*)Titles[ds] << " Tree..."; 
    cout.flush();
    // Print the events in a list
    //TEventList * list = new TEventList(TString("list")+(long)ds,"Zenith[10]>90");
		//list->SetDirectory(0);
		////Lets play with this tree
		//TTreePlayer * tPlay = ((TTree*)Trees[ds])->GetPlayer();
		////Assign selection to tree
		//tPlay->CompileVariables("",Sel.GetTitle());
		////Now run the TreePlayer loop (5 means assign to elist)
		//tPlay->EntryLoop(5,list);
    
		// Change to tree directory
		//gROOT->cd(((TTree*)Trees[ds])->GetDirectory()->GetName());

		// Try dummy drawing before
		((TTree*)Trees[ds])->Project("hdummy","Zenith[1]");
		// Draw into list
		((TTree*)Trees[ds])->Draw((char*)(TString(">> list")+(long)ds),
		                                   Sel.GetTitle(),"");
		//gDirectory->ls(0);
    // Get the list
    TEventList * list = 
		  (TEventList*)gDirectory->Get((char*)(TString("list")+(long)ds));
		if ( ! list ) {
		 cout << "\n NO LIST !!!\n";
		 break;
		}
		list->SetDirectory(0);
    // Assign this list to the Tree
      ((TTree*)Trees[ds])->SetEventList(list);
    if (! Trees[ds]->GetEventList()) {
       cout << "Event list not found !!! \n"; 
       break ;
    } else {
      cout << "\t " << list->GetN() << " Events\n";
    }
  }
}

void CalcPsi( int nentries=0 ){
  // Calculate Psi and the mean for a number of variables. 

	if ( ! Trees[0] ) Init("UPDATE");

  cout << "Loading $ROOTSYS/lib/libPhysics.so ...";
  cout.flush();
  gROOT->LoadMacro("$ROOTSYS/lib/libPhysics.so"); 
  cout << " done\n";
  
  const int nfit = 12;

  int DoFor[nfit] = { kFALSE,	//Linefit
                      kFALSE,	//DirectWalk
                      kFALSE,	//Like1Iter(LF)
                      kTRUE,	//Like1Iter(DW)
                      kFALSE,	//Dipole
                      kFALSE,	//Tensor of inertia
                      kTRUE,	//Like16Iter
                      kFALSE,	//Energy
                      kFALSE,	//Cascade
                      kFALSE,	//Bayesian
                      kTRUE,	//Zylinder-Like16Iter
                      kTRUE,	//DirectWalk II
                     };

  
  // Array to store variables
  float zen[nfit];
  float azi[nfit];

  // Debug histos
  TH1F * hpsi;
  TH1F * hmean;
  TCanvas * myC = new TCanvas("Meiner","Einer");
  myC->Divide(2,1);
  
  //Loop over samles
  for ( int spl = (int)cexp ; spl <= (int)coinc ; spl++ ) {
    cout << "Calculating for " << (char*)Titles[spl] << " data \n";
    // Point tree to array
    Trees[spl]->SetBranchAddress("Zenith",zen);
    Trees[spl]->SetBranchAddress("Azimuth",azi);
    // See if branch is already defined
    TBranch * newBranch = Trees[spl]->GetBranch("Psi");
    if ( newBranch ) {
      cout << "Branch is already defined !\n";
      break ;
    }

    // Debug histo
    hpsi = new TH1F("hpsi","PSI",97,0,3.14159);
    hmean = new TH1F("hpsi","PSI",97,0,180);
    myC->cd(1);
    hpsi->Draw();
    myC->cd(2);
    hmean->Draw();
    
    // We need to store the counted numbers
    struct psi_t { float psi, zen, azi; } mean = { 0,0,0 };
    
    // Now get it new
    newBranch=Trees[spl]->Branch("Mean",&mean,"psi/F:zen/F:azi/F");
    //See how many entries to process - default = all
    int maxentries =0;
    if ( ! nentries )  maxentries = Trees[spl]->GetEntries();
      else maxentries = (Trees[spl]->GetEntries() > nentries ) ?
               nentries : Trees[spl]->GetEntries() ;
    // Used fits
    int ufit =0; for ( int i=0 ; i < nfit; i++) if ( DoFor[i] ) ufit++;
    cout << " Found " << ufit << " fits to use\n";
    //Loop over entries
    for ( long ientry=0 ; ientry < maxentries ; ientry++){
      Trees[spl]->GetEntry(ientry) ;
      //Reset
      mean.zen =0; mean.azi= 0; mean.psi=0;
      // Now loop over fits
      for ( int i=0 ; i < nfit; i++){
        // First get the mean
        if  ( DoFor[i] ) {
          mean.zen += zen[i];
          mean.azi += azi[i];
        }
      }
      // Get the mean vector
      mean.zen /= 1.*ufit; mean.azi /= 1.*ufit;
      TVector3 Mean(0,0,1);
      Mean.SetTheta(mean.zen*3.14156/180);
      Mean.SetPhi(mean.azi*3.14156/180);
      Mean.SetMag(1);
      TVector3 Fit(0,0,1);
      // Now get psi
      for ( int i=0 ; i < nfit; i++){
        if  ( DoFor[i] ) {
          // First get the vector
          Fit.SetTheta(zen[i]*3.14156/180);
          Fit.SetPhi(azi[i]*3.14156/180);
          Fit.SetMag(1);
//          cout << "mean = " << mean.zen << ":" << mean.azi << "\n";
//          cout << "fit = " << zen[i] << ":" << azi [i]<< "\n";
//          cout << "Psi = " << Fit.Angle(Mean)*180/3.14156 << "\n";
          mean.psi += Fit.Angle(Mean);
        }
      }
			mean.psi /= 1.*ufit;
      // fill the debug histo
      hpsi.Fill(mean.psi);
      hmean.Fill(mean.zen);
      newBranch->Fill();
      // be a little verbose
      if ( ientry % 1000 == 0 ) {
        cout << ientry << " Entries \r";
        cout.flush();
        myC->cd(1);
        hpsi->Draw();
        myC->cd(2);
        hmean->Draw();
        myC->Update();
      }
    }
  }
}



void GaryCuts(){
  if ( ! fexp ) loadfile("clean");
  if ( ! ftree("fcoinc") ) loadfriends();

  //Define the cuts

  TCut bayes("Zenith[9]>90");
  TCut ldirb("Ldirb[6]>100");
  TCut smooth("abs(Smootallphit[6])<0.3");
  TCut psi("abs(Zenith[0]-Zenith[6])<20");
  TCut ndirc("Ndirc[6]>7");

  TCanvas * gc = new TCanvas("GCut","Gary's Cuts");
  gc->Divide(2,3);
  gc->cd(1);
  plotall("Zenith[9]","",0,180,90,real);
  gc->cd(2);
  plotall("Ldirb[6]",bayes,0,200,100,real);
  gc->cd(3);
  plotall("Smootallphit[6]",bayes&&ldirb,-1,1,100,real);
  gc->cd(4);
  plotall("Zenith[0]-Zenith[6]",bayes&&ldirb&&smooth,-100,100,100,real);
  gc->cd(5);
  plotall("Ndirc[6]",bayes&&ldirb&&smooth&&psi,0,40,40,real);
  gc->cd(6);
  plotall("Zenith[6]",bayes&&ldirb&&smooth&&psi&&ndirc,0,180,90,real);

}

void HartillCuts(){
  if ( ! fexp ) loadfile("clean");
  if ( ! ftree("fcoinc") ) loadfriends();

  //Define the cuts

  TCut costh("Zenith[6]>95");
  TCut lupdn("Jkchi[6]-Jkprob[6] < 10");
  TCut smooth("abs(Smootallphit[6])<0.3");
  TCut psi("abs(Zenith[0]-Zenith[6])<20");
  TCut ndirc("Ndirc[6]>7");

  TCanvas * gc = new TCanvas("GCut","Gary's Cuts");
  gc->Divide(2,3);
  gc->cd(1);
  plotall("Zenith[9]","",0,180,90,real);
  gc->cd(2);
  plotall("Ldirb",bayes,0,200,100,real);
  gc->cd(3);
  plotall("Smootallphit[6]",bayes&&ldirb,-1,1,100,real);
  gc->cd(4);
  plotall("Zenith[0]-Zenith[6]",bayes&&ldirb&&smooth,-100,100,100,real);
  gc->cd(5);
  plotall("Ndirc[6]",bayes&&ldirb&&smooth&&psi,0,40,40,real);
  gc->cd(6);
  plotall("Zenith[6]",bayes&&ldirb&&smooth&&psi&&ndirc,0,180,90,real);

}