//////////////////////////////////////////////////////////
// This class has been automatically generated on
// Fri Jan 12 17:42:22 2024 by ROOT version 6.30/02
// from TTree events/AMPT tree
// found on file: ampt01.root
//////////////////////////////////////////////////////////

#ifndef Correlator_h
#define Correlator_h

#include <TROOT.h>
#include <TChain.h>
#include <TFile.h>
#include <TComplex.h>
#include <fstream>
// Header file for the classes stored in the TTree if any.
//#include "UEvent.h"
#include "TObject.h"
//#include "UParticle.h"

class Correlator {
public :
   TTree          *fChain;   //!pointer to the analyzed TTree or TChain
   Int_t           fCurrent; //!current Tree number in a TChain

// Fixed size dimensions of array or collections stored in the TTree if any.
   static constexpr Int_t kMaxfParticles = 854;

   // Declaration of leaf types
 //UEvent          *event;
   UInt_t          fUniqueID;
   UInt_t          fBits;
   Int_t           fEventNr;
   Double_t        fB;
   Double_t        fPhi;
   Int_t           fNes;
   Int_t           fStepNr;
   Double_t        fStepT;
   Int_t           fNpa;
   TString         fComment;
   Int_t           fParticles_;
   UInt_t          fParticles_fUniqueID[kMaxfParticles];   //[fParticles_]
   UInt_t          fParticles_fBits[kMaxfParticles];   //[fParticles_]
   Int_t           fParticles_fIndex[kMaxfParticles];   //[fParticles_]
   Int_t           fParticles_fPdg[kMaxfParticles];   //[fParticles_]
   Int_t           fParticles_fStatus[kMaxfParticles];   //[fParticles_]
   Int_t           fParticles_fParent[kMaxfParticles];   //[fParticles_]
   Int_t           fParticles_fParentDecay[kMaxfParticles];   //[fParticles_]
   Int_t           fParticles_fMate[kMaxfParticles];   //[fParticles_]
   Int_t           fParticles_fDecay[kMaxfParticles];   //[fParticles_]
   Int_t           fParticles_fChild[kMaxfParticles][2];   //[fParticles_]
   Double32_t      fParticles_fPx[kMaxfParticles];   //[fParticles_]
   Double32_t      fParticles_fPy[kMaxfParticles];   //[fParticles_]
   Double32_t      fParticles_fPz[kMaxfParticles];   //[fParticles_]
   Double32_t      fParticles_fE[kMaxfParticles];   //[fParticles_]
   Double32_t      fParticles_fX[kMaxfParticles];   //[fParticles_]
   Double32_t      fParticles_fY[kMaxfParticles];   //[fParticles_]
   Double32_t      fParticles_fZ[kMaxfParticles];   //[fParticles_]
   Double32_t      fParticles_fT[kMaxfParticles];   //[fParticles_]
   Double32_t      fParticles_fWeight[kMaxfParticles];   //[fParticles_]

   // List of branches
   TBranch        *b_event_fUniqueID;   //!
   TBranch        *b_event_fBits;   //!
   TBranch        *b_event_fEventNr;   //!
   TBranch        *b_event_fB;   //!
   TBranch        *b_event_fPhi;   //!
   TBranch        *b_event_fNes;   //!
   TBranch        *b_event_fStepNr;   //!
   TBranch        *b_event_fStepT;   //!
   TBranch        *b_event_fNpa;   //!
   TBranch        *b_event_fComment;   //!
   TBranch        *b_event_fParticles_;   //!
   TBranch        *b_fParticles_fUniqueID;   //!
   TBranch        *b_fParticles_fBits;   //!
   TBranch        *b_fParticles_fIndex;   //!
   TBranch        *b_fParticles_fPdg;   //!
   TBranch        *b_fParticles_fStatus;   //!
   TBranch        *b_fParticles_fParent;   //!
   TBranch        *b_fParticles_fParentDecay;   //!
   TBranch        *b_fParticles_fMate;   //!
   TBranch        *b_fParticles_fDecay;   //!
   TBranch        *b_fParticles_fChild;   //!
   TBranch        *b_fParticles_fPx;   //!
   TBranch        *b_fParticles_fPy;   //!
   TBranch        *b_fParticles_fPz;   //!
   TBranch        *b_fParticles_fE;   //!
   TBranch        *b_fParticles_fX;   //!
   TBranch        *b_fParticles_fY;   //!
   TBranch        *b_fParticles_fZ;   //!
   TBranch        *b_fParticles_fT;   //!
   TBranch        *b_fParticles_fWeight;   //!

   Correlator(TTree *tree=0);
   virtual ~Correlator();
   virtual Int_t    Cut(Long64_t entry);
   virtual Int_t    GetEntry(Long64_t entry);
   virtual Long64_t LoadTree(Long64_t entry);
   virtual void     Init(TTree *tree);
   virtual void     Loop();
   virtual Bool_t   Notify();
   virtual void     Show(Long64_t entry = -1);
   
   //User Defined Functions
   
   virtual void    	Cosmetics();
   virtual void    	CalculateQvectors();
   virtual TComplex	Q(Int_t ,Int_t );
   virtual TComplex     Two(Int_t,Int_t);
   
   //User Defined Variables
   
      
 
Bool_t bUseWeights = kTRUE; 

// Cross-check with nested loops:
Bool_t bEvaluateNestedLoops = kTRUE;
   
   // Pick up randomly some harmonics:
//static const Int_t h1=1, h2=3, h3=5, h4=0, h5=-2, h6=-4, h7=-1, h8=-6;
static const Int_t h1=1, h2=2;//, h3=5, h4=0, h5=-2, h6=-4, h7=-1, h8=-6;

// Book Q-vector components: 
//static const Int_t sum = (h1<0?-1*h1:h1)+(h2<0?-1*h2:h2)+(h3<0?-1*h3:h3)+(h4<0?-1*h4:h4)
 //               + (h5<0?-1*h5:h5)+(h6<0?-1*h6:h6)+(h7<0?-1*h7:h7)+(h8<0?-1*h8:h8);
static const Int_t sum = (h1<0?-1*h1:h1)+(h2<0?-1*h2:h2);
//const Int_t maxCorrelator = 2; // We will not go beyond 8-p correlations
//const Int_t maxHarmonic = sum+1; 
//const Int_t maxPower = maxCorrelator+1; 

 static const Int_t maxCorrelator = 2; // We will not go beyond 8-p correlations
 static const Int_t maxHarmonic = sum+1; 
 static const Int_t maxPower = maxCorrelator+1; 
TComplex Qvector[maxHarmonic][maxPower]; // All needed Q-vector components

// Store the final results here:
//  Remark: [2][maxCorrelator] => [Cos,Sin][<2>,<3>,<4>,<5>,<6>,<7>,<8>]
TProfile *correlations[2][maxCorrelator] = {{NULL}}; // Correlations calculated from Q-vector componets using standalone formulas
// Use or not the non-unit particle weights above:

};

#endif

#ifdef Correlator_cxx
Correlator::Correlator(TTree *tree) : fChain(0) 
{
// if parameter tree is not specified (or zero), connect the file
// used to generate this class and read the Tree.
   if (tree == 0) {
      TFile *f = (TFile*)gROOT->GetListOfFiles()->FindObject("ampt01-100-1-1.root");
      if (!f || !f->IsOpen()) {
         f = new TFile("ampt01-100-1-1.root");
      }
      f->GetObject("events",tree);

   }
   Init(tree);
}

Correlator::~Correlator()
{
   if (!fChain) return;
   delete fChain->GetCurrentFile();
}

Int_t Correlator::GetEntry(Long64_t entry)
{
// Read contents of entry.
   if (!fChain) return 0;
   return fChain->GetEntry(entry);
}
Long64_t Correlator::LoadTree(Long64_t entry)
{
// Set the environment to read one entry
   if (!fChain) return -5;
   Long64_t centry = fChain->LoadTree(entry);
   if (centry < 0) return centry;
   if (fChain->GetTreeNumber() != fCurrent) {
      fCurrent = fChain->GetTreeNumber();
      Notify();
   }
   return centry;
}

void Correlator::Init(TTree *tree)
{
   // The Init() function is called when the selector needs to initialize
   // a new tree or chain. Typically here the branch addresses and branch
   // pointers of the tree will be set.
   // It is normally not necessary to make changes to the generated
   // code, but the routine can be extended by the user if needed.
   // Init() will be called many times when running on PROOF
   // (once per file to be processed).

   // Set branch addresses and branch pointers
   if (!tree) return;
   fChain = tree;
   fCurrent = -1;
   fChain->SetMakeClass(1);

   fChain->SetBranchAddress("fUniqueID", &fUniqueID, &b_event_fUniqueID);
   fChain->SetBranchAddress("fBits", &fBits, &b_event_fBits);
   fChain->SetBranchAddress("fEventNr", &fEventNr, &b_event_fEventNr);
   fChain->SetBranchAddress("fB", &fB, &b_event_fB);
   fChain->SetBranchAddress("fPhi", &fPhi, &b_event_fPhi);
   fChain->SetBranchAddress("fNes", &fNes, &b_event_fNes);
   fChain->SetBranchAddress("fStepNr", &fStepNr, &b_event_fStepNr);
   fChain->SetBranchAddress("fStepT", &fStepT, &b_event_fStepT);
   fChain->SetBranchAddress("fNpa", &fNpa, &b_event_fNpa);
   fChain->SetBranchAddress("fComment", &fComment, &b_event_fComment);
   fChain->SetBranchAddress("fParticles", &fParticles_, &b_event_fParticles_);
   fChain->SetBranchAddress("fParticles.fUniqueID", fParticles_fUniqueID, &b_fParticles_fUniqueID);
   fChain->SetBranchAddress("fParticles.fBits", fParticles_fBits, &b_fParticles_fBits);
   fChain->SetBranchAddress("fParticles.fIndex", fParticles_fIndex, &b_fParticles_fIndex);
   fChain->SetBranchAddress("fParticles.fPdg", fParticles_fPdg, &b_fParticles_fPdg);
   fChain->SetBranchAddress("fParticles.fStatus", fParticles_fStatus, &b_fParticles_fStatus);
   fChain->SetBranchAddress("fParticles.fParent", fParticles_fParent, &b_fParticles_fParent);
   fChain->SetBranchAddress("fParticles.fParentDecay", fParticles_fParentDecay, &b_fParticles_fParentDecay);
   fChain->SetBranchAddress("fParticles.fMate", fParticles_fMate, &b_fParticles_fMate);
   fChain->SetBranchAddress("fParticles.fDecay", fParticles_fDecay, &b_fParticles_fDecay);
   fChain->SetBranchAddress("fParticles.fChild[2]", fParticles_fChild, &b_fParticles_fChild);
   fChain->SetBranchAddress("fParticles.fPx", fParticles_fPx, &b_fParticles_fPx);
   fChain->SetBranchAddress("fParticles.fPy", fParticles_fPy, &b_fParticles_fPy);
   fChain->SetBranchAddress("fParticles.fPz", fParticles_fPz, &b_fParticles_fPz);
   fChain->SetBranchAddress("fParticles.fE", fParticles_fE, &b_fParticles_fE);
   fChain->SetBranchAddress("fParticles.fX", fParticles_fX, &b_fParticles_fX);
   fChain->SetBranchAddress("fParticles.fY", fParticles_fY, &b_fParticles_fY);
   fChain->SetBranchAddress("fParticles.fZ", fParticles_fZ, &b_fParticles_fZ);
   fChain->SetBranchAddress("fParticles.fT", fParticles_fT, &b_fParticles_fT);
   fChain->SetBranchAddress("fParticles.fWeight", fParticles_fWeight, &b_fParticles_fWeight);
   Notify();
}

Bool_t Correlator::Notify()
{
   // The Notify() function is called when a new file is opened. This
   // can be either for a new TTree in a TChain or when when a new TTree
   // is started when using PROOF. It is normally not necessary to make changes
   // to the generated code, but the routine can be extended by the
   // user if needed. The return value is currently not used.

   return kTRUE;
}

void Correlator::Show(Long64_t entry)
{
// Print contents of entry.
// If entry is not specified, print current entry
   if (!fChain) return;
   fChain->Show(entry);
}
Int_t Correlator::Cut(Long64_t entry)
{
// This function may be called from Loop.
// returns  1 if entry is accepted.
// returns -1 otherwise.
   return 1;
}

//user Defined Functions

void Correlator::Cosmetics()
{
 // Book everything here.
  
 for(Int_t cs=0;cs<2;cs++) 
 {
  for(Int_t c=0;c<maxCorrelator;c++)
  {
   correlations[cs][c] = new TProfile("","",1,0.,1.);
   correlations[cs][c]->Sumw2();
   
  } // end of for(Int_t c=0;c<maxCorrelator;c++)
 } // end of for(Int_t cs=0;cs<2;cs++) 

} // void Cosmetics()

//==========================================================================================================================

void Correlator::CalculateQvectors()
{
 // Calculate Q-vectors.

 // a) Make sure all Q-vectors are initially zero;
 // b) Calculate Q-vectors for available angles and weights. 

 // a) Make sure all Q-vectors are initially zero:
 for(Int_t h=0;h<maxHarmonic;h++)
 {
  for(Int_t p=0;p<maxPower;p++)
  {
   Qvector[h][p] = TComplex(0.,0.);
  } //  for(Int_t p=0;p<maxPower;p++)
 } // for(Int_t h=0;h<maxHarmonic;h++)

 // b) Calculate Q-vectors for available angles and weights: 
 
 
   if (fChain == 0) return;

   Long64_t nentries = fChain->GetEntriesFast();
   Long64_t nbytes = 0, nb = 0;

//Event Loop
   for (Long64_t jentry=0; jentry<nentries;jentry++) {
      Long64_t ientry = LoadTree(jentry);
      if (ientry < 0) break;
      nb = fChain->GetEntry(jentry);   nbytes += nb;
      
 		for(Int_t i=0;i<fNpa;i++) // loop over particles
 			{
 			
 			if (fParticles_fPy[i]==0 && fParticles_fPx[i] ==0) continue; //exclude spectators
                        Int_t pid1=fParticles_fPdg[i]; 
                        if (abs(pid1)!=2212)continue;
 			Double_t dPhi = TMath::ATan2(fParticles_fPy[i],fParticles_fPx[i]); // particle angle
 			//Double_t mwPhi = 1.; // particle weight
 			Double_t wPhi=1.;
 			Double_t mwPhi = TMath::Sqrt(fParticles_fPy[i]*fParticles_fPy[i]+fParticles_fPx[i]*fParticles_fPx[i]); // particle weight
 			Double_t wPhiToPowerP = 1.; // particle weight raised to power p
 			
 			
  			
  			if(bUseWeights){wPhi=mwPhi ;}   
  				for(Int_t h=0;h<maxHarmonic;h++)
 				{
   	   				for(Int_t p=0;p<maxPower;p++)
   					{
    					if(bUseWeights){wPhiToPowerP = pow(wPhi,p);}
    					
    					Qvector[h][p] += TComplex(wPhiToPowerP*TMath::Cos(h*dPhi),wPhiToPowerP*TMath::Sin(h*dPhi));
    					
    					} //  for(Int_t p=0;p<maxPower;p++)
  				} // for(Int_t h=0;h<maxHarmonic;h++)
  				
  				
 			} //  for(Int_t i=0;i<nParticles;i++) // loop over particles
}

} // void CalculateQvectors()

//==========================================================================================================================

TComplex Correlator::Q(Int_t n, Int_t p)
{
 // Using the fact that Q{-n,p} = Q{n,p}^*. 
 
 if(n>=0){return Qvector[n][p];} 
 return TComplex::Conjugate(Qvector[-n][p]);
 
} // TComplex Q(Int_t n, Int_t p)

//==========================================================================================================================

//==========================================================================================================================


TComplex Correlator::Two(Int_t n1, Int_t n2)
{
 // Generic two-particle correlation <exp[i(n1*phi1+n2*phi2)]>.

 TComplex two = Q(n1,1)*Q(n2,1)-Q(n1+n2,2);

 return two;

} // TComplex Two(Int_t n1, Int_t n2)

//=======================================================================================================================

#endif // #ifdef Correlator_cxx