// 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 "TMath.h"
#include "TClonesArray.h"
#include "/panda/pandaroot/lmd/LmdQA/PndLmdTrackQ.h"
#include <iostream>
#include<algorithm>
using namespace std;
bool myfunction_ch_trio (pair<vector<int>,vector<int> > i,pair<vector<int>,vector<int > > j) { 
  vector<int> datai = i.second;
  vector<int> dataj = j.second;
  return (datai[0]>dataj[0]); 
}
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_samples(TString path="/panda/myResults/DPMoutputDec2013/all_15/results/", double cs_tot = 50.26){
  int simevents = 1e7;
  //    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 lists of all final states particles -------------------------------------
    TString dates[3] = {"18122013","19122013","20122013"}; //15 GeV/c
  //  TString dates[3] = {"25122013","26122013","27122013"}; //1.5 GeV/c
 vector<pair<vector<int>,int> > fchparts[3];

  for(int id=0;id<3;id++){//go over each sample
  TString filename2 = path+"/TOT_Lumi_TrksQA_";//one sample for el+inel mode
  filename2 +=dates[id];
  filename2 +=".root";
  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)
  if(nEvents!=simevents) cout<<"! Attention ! number of events in tree isn't equal number of events asked to be analysed"<<endl;
  //    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<simevents; 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); 
 
    ///c) save channel if new, increase counter if old -------------------------------------
    for(std::vector<pair<std::vector<int>,int> >::iterator itl = fchparts[id].begin(); itl != fchparts[id].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[id].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[id])/sizeof(fchparts[id][0].second);
 Int_t nvalsch =  ch_parts;
 std::sort (fchparts[id].begin(), fchparts[id].end(), myfunction_ch); 
 cout<<"There are "<<nvalsch<<" channels!"<<endl;
 }//end go over each sample
  //[END] 1. ------------------------------------------------------------------------------------------------

  ///2. Make avarage of lists of all final states particles (final states list can be slightly different in each smaple) ----
  std::vector<std::vector<int> > fchpartsGL; //combined list of all possible final states
  std::vector<int> fchpartsGL_vals[3]; //cross-sections values for combined list of all possible final states in each sample
  std::vector<int> fchpartsGL_mean; //cross-sections means for combined list of all possible final states
  std::vector<int> fchpartsGL_sigma; //cross-sections sigmas for combined list of all possible final states
  ///a) make a combined list ----------------------------------------------------------------------------------------------------------
  for(int id1=0;id1<3;id1++){//1st loop over each sample
    for(std::vector<pair<std::vector<int>,int> >::iterator itl = fchparts[id1].begin(); itl != fchparts[id1].end(); ++itl){ //loop over sub-channels in each list
      bool addch = true;
      for(std::vector<std::vector<int> >::iterator itgl = fchpartsGL.begin(); itgl != fchpartsGL.end(); ++itgl){ //loop over already saved sub-channels in global list
	if((*itgl)==(*itl).first) addch = false;
	//	if((*itgl)!=(*itl).first) addch = true;
      }// end loop over already saved sub-channels in global list
      if(addch){ 
	fchpartsGL.push_back((*itl).first);
	fchpartsGL_vals[0].push_back(0);  fchpartsGL_vals[1].push_back(0); 	fchpartsGL_vals[2].push_back(0);
	fchpartsGL_mean.push_back(0);   fchpartsGL_sigma.push_back(0); 
      }
    }//end loop over sub-channels in each list
  }//end 1nd loop over each sample

 int glcount = 0;
  for(std::vector<std::vector<int> >::iterator itgl = fchpartsGL.begin(); itgl != fchpartsGL.end(); ++itgl){ //loop over saved sub-channels in global list
    for(int id1=0;id1<3;id1++){//1st loop over each sample
       for(std::vector<pair<std::vector<int>,int> >::iterator itl = fchparts[id1].begin(); itl != fchparts[id1].end(); ++itl){ //loop over sub-channels in each list
	 if((*itgl)==(*itl).first){
	     fchpartsGL_vals[id1][glcount] = (*itl).second;
	     //     continue;
	   }
       }//end loop over sub-channels in each list
       // if(glcount==1) cout<<" at sample "<<id1<<" for ch#"<<glcount<<" there are "<<fchpartsGL_vals[id1][glcount]<<" events"<<endl;
    }//end 1nd loop over each sample
       glcount++;
  }// end loop over saved sub-channels in global list

  ///b) calculate avarage means of combined list ----------------------------------------------------
  for(int iglc = 0;iglc<glcount;iglc++){//loop over saved sub-channels in global list
    for(int id1=0;id1<3;id1++){//1st loop over each sample
      fchpartsGL_mean[iglc]  += fchpartsGL_vals[id1][iglc];
    }//end 1st loop over each sample
    fchpartsGL_mean[iglc] *=1./3;
  }//endl loop over saved sub-channels in global list

  ///c) calculate avarage sigmas of combined list ----------------------------------------------------
  for(int iglc = 0;iglc<glcount;iglc++){//loop over saved sub-channels in global list
    for(int id1=0;id1<3;id1++){//1st loop over each sample
      fchpartsGL_sigma[iglc]  +=(fchpartsGL_vals[id1][iglc]-fchpartsGL_mean[iglc])*(fchpartsGL_vals[id1][iglc]-fchpartsGL_mean[iglc]);
    }//end 1st loop over each sample
    fchpartsGL_sigma[iglc] *=1./3;
    //    double  sigma = TMath::St(fchpartsGL_sigma[iglc]);
    fchpartsGL_sigma[iglc] = sqrt(fchpartsGL_sigma[iglc]);
  }//endl loop over saved sub-channels in global list

 //[END] 2. ------------------------------------------------------------------------------------------------

  /// 3. make new combined list: with channels and contributions info ------------------------
  std::vector<pair<std::vector<int>,std::vector<int> > > fchinfo;
  int glc = 0;
 for(std::vector<std::vector<int> >::iterator itl = fchpartsGL.begin(); itl != fchpartsGL.end(); ++itl){ //loop over saved sub-channels
   vector<int> vcs;
   vcs.push_back(fchpartsGL_mean[glc]);
   vcs.push_back(fchpartsGL_sigma[glc]);
   vector<int> vch = (*itl);
   fchinfo.push_back(make_pair(vch,vcs));
   glc++;
 }// end loop over saved sub-channels
 std::sort(fchinfo.begin(),fchinfo.end(), myfunction_ch_trio); 
  ///---------------------------------------------------------------------------------------------------------

  ///4 print avaraged results ----------------------------------------------------------------------------
 double tot_sum_ck=0;
 double tot_sum_ck_err=0;
 glc = 0;
 for(std::vector<pair<std::vector<int>,std::vector<int> > >::iterator ich = fchinfo.begin(); ich != fchinfo.end(); ++ich){// loop over saved & ranged sub-channels
   vector<int> subch = (*ich).first;
   vector<int> data = (*ich).second;
   double cs_DPM_mn = double(data[0])*cs_tot/simevents;
   double cs_DPM_sig = double(data[1])*cs_tot/simevents;
   tot_sum_ck +=cs_DPM_mn;
    tot_sum_ck_err +=cs_DPM_sig;
  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);
      }
    }
    //double cs_DPM_mn = fchpartsGL_mean[glcount]*cs_tot/simevents;
    //    double cs_DPM_sig = fchpartsGL_sigma[glcount]*cs_tot/simevents;

    if(cs_DPM_mn*1e3<1e-1){
      cout<<" ("<<cs_DPM_mn*1e6<<" +/- "<<cs_DPM_sig*1e6<<" nb): ";
    }
    else{
      if(cs_DPM_mn*1e3>1e3){
    	cout<<" ("<<cs_DPM_mn<<" +/- "<<cs_DPM_sig<<" mb): ";
      }
      else{
    	cout<<" ("<<cs_DPM_mn*1e3<<" +/- "<<cs_DPM_sig*1e3<<" mub): ";
      }
    }
    //  cout<<" with "<<fchpartsGL_mean[glcount]<<" events ("<<fchpartsGL_vals[0][glcount]<<", "<<fchpartsGL_vals[1][glcount]<<", "<<fchpartsGL_vals[2][glcount]<<")"<<endl;
    cout<<" "<<endl;
   glc++;
 }//end loop over saved & ranged sub-channels
 cout<<" TOT sum for "<<glcount<<" channels: "<<tot_sum_ck<<" +/- "<<tot_sum_ck_err<<" mb"<<endl;
  ///----------------------------------------------------------------------------------------------------------


  // ///3 print avaraged results ----------------------------------------------------------------------------
  // double tot_sum_ck=0;
  // double tot_sum_ck_err=0;
  // glcount = 0;
  // //  for(std::vector<pair<std::vector<int>,int> >::iterator itl = fchparts[id].begin(); itl != fchparts[id].end(); ++itl){ //loop over saved sub-channels
  // for(std::vector<std::vector<int> >::iterator itl = fchpartsGL.begin(); itl != fchpartsGL.end(); ++itl){ //loop over saved sub-channels
  //   //    vector<int> subch = (*itl).first;
  //   //    if(fchpartsGL_mean[glcount]<2) continue;
  //   vector<int> subch = (*itl);
  //   double cs_DPM_mn = double(fchpartsGL_mean[glcount])*cs_tot/simevents;
  //   double cs_DPM_sig = double(fchpartsGL_sigma[glcount])*cs_tot/simevents;
  //   tot_sum_ck +=cs_DPM_mn;
  //   tot_sum_ck_err +=cs_DPM_sig;
  //   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);
  //     }
  //   }
  //   //double cs_DPM_mn = fchpartsGL_mean[glcount]*cs_tot/simevents;
  //   //    double cs_DPM_sig = fchpartsGL_sigma[glcount]*cs_tot/simevents;

  //   if(cs_DPM_mn*1e3<1e-1){
  //     cout<<" ("<<cs_DPM_mn*1e6<<" +/- "<<cs_DPM_sig*1e6<<" nb): ";
  //   }
  //   else{
  //     if(cs_DPM_mn*1e3>1e3){
  //   	cout<<" ("<<cs_DPM_mn<<" +/- "<<cs_DPM_sig<<" mb): ";
  //     }
  //     else{
  //   	cout<<" ("<<cs_DPM_mn*1e3<<" +/- "<<cs_DPM_sig*1e3<<" mub): ";
  //     }
  //   }
  //   //  cout<<" with "<<fchpartsGL_mean[glcount]<<" events ("<<fchpartsGL_vals[0][glcount]<<", "<<fchpartsGL_vals[1][glcount]<<", "<<fchpartsGL_vals[2][glcount]<<")"<<endl;
  //   cout<<" "<<endl;
  //   glcount++;
  // }//end loop over saved sub-channels
  // cout<<" TOT sum for "<<glcount<<" channels: "<<tot_sum_ck<<" +/- "<<tot_sum_ck_err<<" mb"<<endl;
  // ///----------------------------------------------------------------------------------------------


  //  }//end go over each sample
}