// macro for analys bkg channels in DPM
// takes TrkQA from LMD files and calculates cross-section for each channel.
// Attention! Should take into account only primary final state!
// author: a.karavdina
// date: 21 January 2014

#include "TString.h"
#include "TFile.h"
#include "TTree.h"
#include "TBranch.h"
#include "TChain.h"
#include "TClonesArray.h"
#include "/panda/pandaroot/lmd/LmdQA/PndLmdTrackQ.h"
#include <iostream>
#include<algorithm>
using namespace std;
bool myfunction_ch (pair<vector<int>,int> i,pair<vector<int>,int> j) { return (i.second>j.second); }
bool myfunction_el (int i,int j) { return (i<j); }
void DPMinel(TString path="/panda/myResults/DPMoutputDec2013/all_15/results/", double cs_tot = 50.26){
  int simevents = 3e7;
  //    const double cs_tot = 50.26;// mb, el+inel @ 15 GeV/c
  // const double cs_tot = 124.63;// mb, el+inel @ 1.5 GeV/c
  ///1. Make list of all final states particles -------------------------------------
  TString filename2 = path+"/TOT_Lumi_TrksQA_el_inel.root";//sum samples for el+inel mode
  //    TString filename2 = path;//one sample for el+inel mode
  TChain tTrkRec("cbmsim");
  TString in=filename2;
  tTrkRec.Add(in);
  //--- Tracks -------------------------------------------------------------------
  TClonesArray* qtrk=new TClonesArray("PndLmdTrackQ");
  tTrkRec.SetBranchAddress("LMDTrackQ",&qtrk);  //Tracks
  //----------------------------------------------------------------------------------
  const int nEvents= tTrkRec.GetEntries();//~(number of simulated events) * (number of tracks per event)
  //    const int nEvents=1e4;//number of simulated events

  /// a) let's make lists of each sub-channel ---------------------------------------------
  vector<pair<vector<int>,int> > fchparts;
  int ch_parts=0;
  for (Int_t j=0; j<nEvents; j++){ // loop over all events
    bool nxtch = true;
    vector<int> chparts;
    int el_parts = 0;
    tTrkRec.GetEntry(j);
    const int nTrksQ = qtrk->GetEntriesFast();//Number of trks/ev
    if(nTrksQ<1) continue;
    for(int iq=0;iq<nTrksQ;iq++){ //loop over all trks in event
      PndLmdTrackQ *trkcur = (PndLmdTrackQ *)qtrk->At(iq);
      int issec = trkcur->GetSecondary();
      if(issec>=0) continue;//primaries only
      int PDGcode = trkcur->GetPDGcode();
      chparts.push_back(PDGcode);
      el_parts++;
    }// end loop over all trks in event
    ///----------------------------------------------------------------------------------------------
    /// b)range elements in the list------------------------------------------------------------
    std::sort (chparts.begin(), chparts.end(), myfunction_el); 
    //    Int_t nvals = sizeof(chparts)/sizeof(chparts[0]); //doesn't work =(
    // Int_t nvals =el_parts;
    // for(int itt=0; itt<nvals*nvals; itt++){//repeat ranging 
    //   int prev_el = chparts[0];
    //   int nc=1;
    //   for(std::vector<int>::iterator isch = chparts.begin()+1; isch != chparts.end(); ++isch){ //loop over particles sub-channel
    // 	int next_el = (*isch);
    // 	if(prev_el>next_el){
    // 	  chparts[nc-1] =next_el;
    // 	  chparts[nc] =prev_el;
    // 	}
    // 	prev_el = next_el;
    // 	nc++;
    //   }// end loop over particles sub-channel
    // }// end repeat ranging 
    ///----------------------------------------------------------------------------------------------
    ///c) save channel if new, increase counter if old -------------------------------------
    for(std::vector<pair<std::vector<int>,int> >::iterator itl = fchparts.begin(); itl != fchparts.end(); ++itl){ //loop over already saved sub-channels
      if(chparts==(*itl).first){
	(*itl).second++;
	nxtch = false;
      }
    }//end loop over already saved sub-channels
    if(nxtch){
      fchparts.push_back(std::make_pair(chparts,1));
      ch_parts++;
    }
  }//end loop over all events
  ///----------------------------------------------------------------------------------------------
///d) range elements in channels vector according to number of events: hights first ----
//  Int_t nvalsch = sizeof(fchparts)/sizeof(fchparts[0].second);
 Int_t nvalsch =  ch_parts;
 cout<<"There are "<<nvalsch<<" channels!"<<endl;

 std::sort (fchparts.begin(), fchparts.end(), myfunction_ch); 

 // for(int itt=0; itt<nvalsch*nvalsch; itt++){//repeat ranging channels
 //   // pair<std::vector<int>,int>  prev_el = fchparts[0];
 //   int nc=1;
 //   int ncitl=0;
 //   for(std::vector<pair<std::vector<int>,int> >::iterator itl = fchparts.begin(); itl != fchparts.end(); ++itl){ //loop over saved sub-channels
 //     ncitl++;
 //     for(std::vector<pair<std::vector<int>,int> >::iterator itl2 = fchparts.begin(); itl2 != (fchparts.end()-ncitl); ++itl2){ //second loop over saved sub-channels
 //       pair<std::vector<int>,int>  prev_el =  (*itl);
 //       pair<std::vector<int>,int>  next_el = (*(itl2+1)); std::next(mylist.begin(),5)
 //     if(prev_el.second<next_el.second){
 //       //       cout<<"exchange of channels with "<<prev_el.second<<" and "<<next_el.second<<" now "<<next_el.second<<" is 1st!"<<endl;
 //       fchparts[nc-1] =next_el;
 //       fchparts[nc] =prev_el;
 //     }
 //     //     prev_el = next_el;
 //     nc++;
 //     }// end second loop over saved sub-channels
 //   }//end loop over saved sub-channels
 // }// end repeat ranging channels
  ///-----------------------------------------------------------------------------------------------------
  ///f) print result ----------------------------------------------------------------------------
  double tot_sum_ck=0;
  for(std::vector<pair<std::vector<int>,int> >::iterator itl = fchparts.begin(); itl != fchparts.end(); ++itl){ //loop over saved sub-channels
    vector<int> subch = (*itl).first;
    tot_sum_ck +=((*itl).second*cs_tot)*1e3/simevents;
  
    for(std::vector<int>::iterator isch = subch.begin(); isch != subch.end(); ++isch){ //loop over particles sub-channel
      switch ((*isch)){//assign human readable values to PDG code
      case 11:
	cout << " e^{-} ";
	break;
      case -11:
	cout << " e^{+} ";
	break;
      case 13:
	cout << " mu^{-} ";
	break;
      case -13:
	cout << " mu^{+} ";
	break;
      case 22:
	cout << " gamma ";
	break;
      case 111:
	cout << " pi^{0} ";
	break;
      case 130:
	cout << " K^{0}_{L} ";
	break;
      case -211:
	cout << " pi^{-} ";
	break;
      case 211:
	cout << " pi^{+} ";
	break;
      case 311:
	cout << " K^{0}";
	break;
      case 310:
	cout << " K^{0}_{S} ";
	break;
      case 321:
	cout << " K^{+} ";
	break;
      case -321:
	cout << " K^{-} ";
	break;
      case -2112:
	cout << " bar{n} ";
	break;
      case 2112:
	cout << " n";
	break;
      case 2212:
	cout << " p ";
	break;
      case -2212:
	cout << " bar{p} ";
	break;
      case 3112:
	cout << " Sigma^{-} ";
	break;
      case 3122:
	cout << " Lambda ";
	break;
      case -3122:
	cout << " bar{Lambda} ";
	break;
      case 3212:
	cout << " Sigma^{0} ";
	break;
      case 3222:
	cout << " Sigma^{+} ";
	break;
      default:
	cout<<" "<<(*isch);
      }
    }
    if(((*itl).second*cs_tot)*1e3/simevents<1e-1){
      cout<<" ("<<((*itl).second*cs_tot)*1e6/simevents<<" nb): ";
    }
    else{
      if(((*itl).second*cs_tot)*1e3/simevents>1e3){
	cout<<" ("<<((*itl).second*cs_tot)/simevents<<" mb): ";
      }
      else{
	cout<<" ("<<((*itl).second*cs_tot)*1e3/simevents<<" mub): ";
      }
    }
    cout<<" "<<endl;
  }//end loop over saved sub-channels
  cout<<" TOT sum: "<<tot_sum_ck*1e-3<<" mb"<<endl;
  ///----------------------------------------------------------------------------------------------
  //[END] 1. ------------------------------------------------------------------------------------------------
}