//************************************************************************//
//
//
//
//
//  To calculate the mean pT distribution for mix events
//
// This calculates the mean pt and the mixed mean pt
// Only for a single Centrality class
//
//
//
//*************************************************************************//


#define PtCorr_mix_cxx
#include "PtCorr_mix.h"
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>


const Int_t nCentbins = 9;



void PtCorr_mix::Loop()
{
//   In a ROOT session, you can do:
//      Root > .L PtCorr_mix.C
//      Root > PtCorr_mix t
//      Root > t.GetEntry(12); // Fill t data members with entry number 12
//      Root > t.Show();       // Show values of entry 12
//      Root > t.Show(16);     // Read and show values of entry 16
//      Root > t.Loop();       // Loop on all entries
//

//     This is the loop skeleton where:
//    jentry is the global entry number in the chain
//    ientry is the entry number in the current Tree
//  Note that the argument to GetEntry must be:
//    jentry for TChain::GetEntry
//    ientry for TTree::GetEntry and TBranch::GetEntry
//
//       To read only selected branches, Insert statements like:
// METHOD1:
//    fChain->SetBranchStatus("*",0);  // disable all branches
//    fChain->SetBranchStatus("branchname",1);  // activate branchname
// METHOD2: replace line
//    fChain->GetEntry(jentry);       //read all branches
//by  b_branchname->GetEntry(ientry); //read only this branch

if (fChain == 0) return;

  // Long64_t nentries = fChain->GetEntriesFast();  //kBigNumber
   Int_t nEvent = fChain->GetEntries();
   //cout << "nentries = " << nentries << endl;
   cout << "nEvent = " << nEvent << endl;
   

   fpw = new TFile("meanptmix.root","RECREATE");
   tree = new TTree("tree","Mean pt plots for mix events");

   tree->Branch("tmptmix",&tmptmix,"tmptmix/F");

   TH1F *hmpt = new TH1F("hmpt"," <<p_{t}>>_{real} ", 200, 0., 2.);
   TH1F *hmpt_mix = new TH1F("hmpt_mix"," <<p_{t}>>_{mix} ", 200, 0., 2.);
   

  // FxtMult Nch cuts

  //Centrality classes      5%    10%  20%  30% 40% 50% 60% 70% 80% 


   //Int_t Nch_cut_up[9]    =  {195,  141, 119, 85, 60, 41, 26, 16, 8 };

   //Int_t Nch_cut_below[9] =  {142,  120, 86,  61, 42, 27, 17, 9,  5 };


   Int_t cent;
   Int_t nevt = 0;  // ---------a counter for that centrality!!------

   const int numevt = 900500;
   //const int numtrk = 400;

   //const int numevt = 10000;
   const int numtrk = 500;

   Int_t totpart[numevt];
   Float_t mpt[numevt];   

   Float_t ptmix[numevt][numtrk];


   Int_t statusnevt[numevt];
   //Int_t statuspt[numevt][numtrk];
 
   for(Int_t i = 0; i < numevt; i++)
     {
       totpart[i] = 0;
       mpt[i] = 0.;
       for(Int_t j = 0; j < numtrk; j++)
	 {
	   ptmix[i][j] = 0.;
	 }
       
     }
  
    //Event Loop.......
   for (Long64_t jentry=0; jentry<nEvent;jentry++) 
     {
       
       fChain->GetEntry(jentry);
       Int_t Nch = fRefMult;
       Double_t pt_sum=0.0;
       
       cout << " Processing event# " << jentry << endl;

       //cout << "Tracks in Event = " << fNtrack << endl;
       
       Int_t counttracks = 0;
    
       if((Nch >= 142) && (Nch < 195) )cent = 0;
        
       else cent =1;

        
       if (cent!=0) continue; // Chooses the 0-5% Centrality class
    
       for(Int_t iTrack = 0; iTrack < fNtrack+1; iTrack++){   // Track Loop

           Double_t px = fPx[iTrack];
	         Double_t py = fPy[iTrack];
	         Double_t pz = fPz[iTrack];
           Double_t p0 = sqrt(px*px+py*py+pz*pz);
           Double_t eta = 0.5*log((p0+pz)/(p0-pz));
	         Double_t pt  = sqrt(px*px + py*py);
           Double_t chg = fCharge[iTrack];
        
            //cout <<px<<py<<pz<< endl;

       
           if(eta>0.0 || eta < -2.0) continue;           // 2 units of eta window
           if( chg==0) continue;                // Only charged particles
           if(pt < 0.15 || pt > 2.0) continue;  // pT  : 0.15 to 2.0 GeV/c
           ptmix[nevt][counttracks] = pt;
           pt_sum += pt;
           counttracks++;
        
        }       //Track Loop End
    

       totpart[nevt] = counttracks;

       cout << "Charged tracks = " << counttracks << endl;

       mpt[nevt] = pt_sum/counttracks;

       hmpt->Fill(mpt[nevt]);

       nevt++;
	//if(nevt>900000) break;

   }//Event Loop ends

/*for(int i = 0; i < numevt; i++){  

    for(int j = 0; j < numtrk; j++){

	     statuspt[i][j] = 0;

	 }

}*/


//----------------------------Mixed Events Analysis----------------------------------
   
   cout << "starting Mixing!!" << endl;
   cout << "Centrality counter = " << nevt << endl;
   Float_t mpt_mix = 0.;

   for(Int_t k = 0; k < nevt; k++){
       
      
       if(k == nevt-2) break;
       
       //Initialization of Event status
       for(Int_t ii = 0; ii < numevt; ii++){
	         
           statusnevt[ii] = 0;
	     
       }

       Float_t pt_mix = 0.;
       Int_t countmix = 0;
       
     lok:
       Int_t event = (gRandom->Rndm())*nevt;
       
       //if(statusnevt[event] == 0)
	 {
	   Int_t part =  (gRandom->Rndm())*totpart[event];
	   
	   //if(statuspt[event][part] == 0){
        if(statusnevt[event] == 0)
	       pt_mix += ptmix[event][part];
	       //statuspt[event][part] =  1;
	       statusnevt[event] = 1;
	       countmix++;

	     }
	 }
     //cout << "countmix = " << countmix << "totpart = " << totpart[k] << endl;
	 if(countmix != totpart[k]) goto lok;

	 //if(TMath::Abs(countmix-totpart[k])<2) goto lok;
       
       
       //cout << " Mixed Event Number ===> " << k << " Completed " << " Tot particles " << countmix << " Tot Particle in actual event " << totpart[k] << endl;

       mpt_mix = pt_mix/countmix;
       tmptmix = mpt_mix;
       hmpt_mix->Fill(mpt_mix);

       //cout <<" mpt mix " << mpt_mix << " mpt real " << mpt[k] << endl;
       


       tree->Fill();

     }// Mixed event loop ends

   hmpt->Draw();
   hmpt_mix->Draw("same");


   fpw->cd();
   hmpt_mix->Write();
   hmpt->Write();
   tree->Write();
   fpw->Close();

}