// macro for analys bkg channels via sum of PDG id
// and print out results
// Idea: estimate the worst case scenario by total cross-section from DPM for all contrinuting to LMD bkg channels
// author: a.karavdina
// date: 23 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 sumPDGcodesALL(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
  }
  //  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
   
      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: "<<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;
  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
   	//now let's make lists of each sub-channel 
	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
	  //   cout<<"read event info with"<<nTrksQ <<" trks"<<endl;
	  if(nTrksQ<1) continue;
	  PndLmdTrackQ *trkcur0 = (PndLmdTrackQ *)qtrk->At(0);
	  int sumPDGev = trkcur0->GetSumEvPDG();
	  if(sumPDGev==(*it).first){//found needed event
	    if(fworstc){ //sum all events from MC
	      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
	    }
	    else{ //check if this event contains reconstructed trks
	      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();
		int trkSt = trkcur->GetTrkRecStatus();
	    	if(trkSt==0) runagain=true; //check only channels, which really contain interesting 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
	    } //end check if this event contains reconstructed trks
	   

	    //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 all events
  }// end loop over channels obtained in previous step

	//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. ------------------------------------------------------------------------------------------------
}