// macro for analys bkg channels via sum of PDG id
// and print out results
// Attention! sum of PDG id corresponds to primary final state only!
// author: a.karavdina
// date: 10 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_el (int i,int j) { return (i<j); }
bool myfunction_ch (pair<vector<int>,int> i,pair<vector<int>,int> j) { return (i.second>j.second); }
void sumPDGcodes(int gl_id=-2212,double fcs_tot=0, TString path="/panda/myResults/DPMoutputDec2013/all_15/results/", bool fworstc = true){
  if(fworstc) cout<<"Count cross-section for the worst case scenario!"<<endl;
  double cs_tot;
  if(fcs_tot!=0){
    cs_tot = fcs_tot;
  }
  else{
    //    cs_tot = 50.26;// mb, el+inel @ 15 GeV/c
    cs_tot = 40.50;// mb, inel @ 15 GeV/c
  }
  //  const double cs_tot = 124.63;// mb, el+inel @ 1.5 GeV/c
  //1. Make list of channels ------------------------------------------------------------------------
  // for particle with PDGid = gl_id
  // -> need loop over event in the tree with "id" and "sumid" branches 
 TString filename1 = path+"/TOT_Lumi_BkgSig_el_inel.root";
  TFile *fin1 = new TFile(filename1,"READ");
  TTree *tbkg1 = (TTree *)fin1->Get("tBkg"); //Bkg data 
  int f_id, f_sumid;
  TBranch *bid = tbkg1->GetBranch("id");
  bid->SetAddress(&f_id);
  TBranch *bsumid = tbkg1->GetBranch("sumid");
  bsumid->SetAddress(&f_sumid);
   int f_issec;
  TBranch *bissec = tbkg1->GetBranch("issecond");
  bissec->SetAddress(&f_issec);

  Int_t nevent = tbkg1->GetEntries();
  Int_t nevent_pri = tbkg1->GetEntries("issecond<0");//only primary events
  vector<pair<int,int> > v_sum_id; //vector of channels: first=id,second=counter
  for (Int_t i=0;i<nevent;i++){
    tbkg1->GetEvent(i);
    if(f_issec>=0) continue;//primaries only
    if(f_id==gl_id){
      //   cout<<" sum_id = "<<f_sumid<<endl;
      //   v_sum_id.push_back(f_sumid);
      bool addit=true;
      for (std::vector<pair<int,int> >::iterator it = v_sum_id.begin(); it != v_sum_id.end(); ++it){
	if((*it).first==f_sumid){
	  addit=false;// add only channels not listed before;
	  (*it).second++;
	}
      }
      if(addit) v_sum_id.push_back(std::make_pair(f_sumid,1));
    }
  }
  std::cout<<"List of channels for particle with ID="<<gl_id<<endl;
 std::sort (v_sum_id.begin(), v_sum_id.end(), myfunction_ch); 
  double contrSum = 0;
  for (std::vector<pair<int,int> >::iterator it = v_sum_id.begin(); it != v_sum_id.end(); ++it){
    contrSum +=100*((*it).second/double(nevent_pri));
    std::cout << ' ' << (*it).first <<" with contribution "<< 100*((*it).second/double(nevent_pri))<<" %  from tot.prim.inel.bkg";
    std::cout << '\n';
  }
  cout<<"TOTAL: "<<contrSum<<endl;
  fin1->Close();
  //[END] 1. ------------------------------------------------------------------------------------------------


  //2. Make list of final state particles for each channel -------------------------------------
  // -> need loop over events in full tree over "LMDTrackQ.fPDGcode" and "LMDTrackQ.fsumID" branches
  TString filename2 = path+"/TOT_Lumi_TrksQA_el_inel.root";//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
  //const int nEvents=1e4;//number of simulated events
    vector<pair<vector<int>,int> > fchparts;
 
   	//now let's make lists of each sub-channel 
    //    for (Int_t j=0; j<nEvents; j++){ // loop over all events
    for (Int_t j=0; j<1000; j++){ // TEST loop over events
      tTrkRec.GetEntry(j);
      const int nTrksQ = qtrk->GetEntriesFast();//Number of trks/ev
      if(nTrksQ<1) continue;
      PndLmdTrackQ *trkcur0 = (PndLmdTrackQ *)qtrk->At(0);
      int sumPDGev = trkcur0->GetSumEvPDG();

      for (std::vector<pair<int,int> >::iterator it = v_sum_id.begin(); it != v_sum_id.end(); ++it){//loop over channels obtained in previous step
   
	  bool nxtch = true;
	  vector<int> chparts;
	  int el_parts=0;
	  if(sumPDGev==(*it).first){//found needed event
	  cout<<" "<<endl;
	    bool runagain = false;
	    for(int iq=0;iq<nTrksQ;iq++){ 
	      PndLmdTrackQ *trkcur = (PndLmdTrackQ *)qtrk->At(iq);
	      int issec = trkcur->GetSecondary();
	      if(issec>=0) continue;//primaries only
	      int PDGcode = trkcur->GetPDGcode();
	      cout<<"PDGcode = "<<PDGcode<<endl;
	      int trkSt = trkcur->GetTrkRecStatus();
	      //   
	      if(fworstc){
		if(PDGcode==gl_id) runagain=true; //check only channels, which really contain interesting particle, doesn't matter was it reconstructed or not -> "worst case"
	      }
	      else{
		if(PDGcode==gl_id && trkSt==0) runagain=true; //check only channels, which really contain interesting & reconstructed particle 
	      }
	    }
	    if(!runagain) continue;
	    for(int iq=0;iq<nTrksQ;iq++){ //loop over all trks in event
	      //     bool addit=true;
	      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

	    //range elements in the list
	    std::sort (chparts.begin(), chparts.end(), myfunction_el); 
	    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));
	  }
    }// end loop over channels obtained in previous step
    }//end loop over all events


    cout<<"ranging elements"<<endl;
	//range elements in channels list
	std::sort (fchparts.begin(), fchparts.end(), myfunction_ch); 
	double tot_sum_ck=0;
	double tot_cs = 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 +=100.*double((*itl).second)/double(nevent_pri);
	  //	  cout<<"List of particles [sumIDs="<<(*it).first<<"]: "<<" ("<<100.*double((*itl).second)/double(nevent_pri)<<" %  from tot.inel.bkg): ";
	  cout<<" ("<<100.*double((*itl).second)/double(nevent_pri)<<" %  from prim.inel.bkg): ";
	  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 3212:
	      cout << " Sigma^{0} ";
	      break;
	    case 3222:
	      cout << " Sigma^{+} ";
	      break;
	    default:
	      cout<<" "<<(*isch);
	    }
	  }
	  cout<<" eff.cs= "<<((*itl).second*cs_tot)/3e1<<" nb"<<endl;
	  tot_cs += ((*itl).second*cs_tot)/3e1;
	}//end loop over saved sub-channels
	if(fworstc)
	  cout<<"The worts case scenario: "<<endl;
	else
	  cout<<"For reconstructed bkg only: "<<endl;
	cout<<"TOTAL eff.cs = "<<tot_cs<<endl;
  //[END] 2. ------------------------------------------------------------------------------------------------
}