[u/mrichter/AliRoot.git] / PWGCF / EBYE / macros / runBalanceFunctionOnHydro.C

//========================================================//
//Macro that reads the output of the hydro calculations 
//(Therminator) using the format that can be found at
//arXiv:1102.0273
//Author: Panos.Christakoglou@nikhef.nl
//========================================================//

//========================================================//
//Event structure
struct StructEvent {
  UInt_t eventID;//unique identifier of the event (random number)
  UInt_t entries;//number of particles of each event
  UInt_t entriesprev;//set to 0 by default
};
//========================================================//

//========================================================//
//Particle structure
//Primordial particles: fathereid==-1 ==> 
//pid==fatherpid==rootpid
struct StructParticle {
  Float_t mass;//the mass of the particle (in GeV)
  Float_t t;//the time coordinate of the particle in fm/c
  Float_t x;//the spacial coordinate x in fm/c
  Float_t y;//the spacial coordinate y in fm/c
  Float_t z;//the spacial coordinate z in fm/c
  Float_t e;//the energy in GeV
  Float_t px;//the x-coordinate of the particle's momentum in GeV
  Float_t py;//the y-coordinate of the particle's momentum in GeV
  Float_t pz;//the z-coordinate of the particle's momentum in GeV
  Int_t decayed;//decay flag (1: decayed, 0: no decay)
  Int_t pid;//PDG of particle
  Int_t fatherpid;//PDG of parent
  Int_t rootpid;//root (primordial) particle PDG number
  Int_t eid;//sequence number in the event
  Int_t fathereid;//parent sequence number in the event
  UInt_t eventid;//unique identifier of the event (random number)
};
//========================================================//

//========================================================//
//Balance function analysis variables
Bool_t gRunShuffling=kFALSE;
Bool_t gRunMixing=kTRUE;
Bool_t gRunMixingWithEventPlane=kFALSE;

Double_t gEtaMin = -0.8;
Double_t gEtaMax = 0.8;
Double_t gPtMin = 0.2;
Double_t gPtMax = 20.0;
//========================================================//


void runBalanceFunctionOnHydro(TString aEventDir = ".", 
			       Int_t aEventFiles = 9) {
  //Macro that reads the themrinator events
  //Author: Panos.Christakoglou@nikhef.nl
  // Time:
  TStopwatch timer;
  timer.Start();

  //========================================================//
  //Load the aliroot libraries
  gSystem->Load("libSTEERBase");
  gSystem->Load("libESD");
  gSystem->Load("libAOD");
  gSystem->Load("libANALYSIS");
  gSystem->Load("libANALYSISalice");
  gSystem->Load("libEventMixing");
  gSystem->Load("libCORRFW");
  gSystem->Load("libPWGTools");
  gSystem->Load("libPWGCFebye");
  //========================================================//

  //========================================================//
  //Configure the bf objects
  AliBalancePsi *bf = new AliBalancePsi();
  bf->SetAnalysisLevel("MC");
  bf->SetShuffle(gRunShuffling);
  bf->SetEventClass("EventPlane");
  bf->SetDeltaEtaMax(TMath::Abs(gEtaMax-gEtaMin));
  bf->InitHistograms();

  AliBalancePsi *bfs = 0x0;
  if(gRunShuffling) {
    bfs = new AliBalancePsi();
    bfs->SetAnalysisLevel("MC");
    bfs->SetEventClass("EventPlane");
    bfs->SetDeltaEtaMax(TMath::Abs(gEtaMax-gEtaMin));
    bfs->InitHistograms();
  }

  AliBalancePsi *bfm = 0x0;
  if(gRunMixing) {
    bfm = new AliBalancePsi();
    bfm->SetAnalysisLevel("MC");
    bfm->SetShuffle(gRunShuffling);
    bfm->SetEventClass("EventPlane");
    bfm->SetDeltaEtaMax(TMath::Abs(gEtaMax-gEtaMin));
    bfm->InitHistograms();
  }
  //========================================================//

  //========================================================//
  //Output objects
  //QA list
  fList = new TList();
  fList->SetName("listQA");
  fList->SetOwner();

  //Balance Function list
  TList *fListBF = new TList();
  fListBF->SetName("listBF");
  fListBF->SetOwner();
  fListBF->Add(bf->GetHistNp());
  fListBF->Add(bf->GetHistNn());
  fListBF->Add(bf->GetHistNpn());
  fListBF->Add(bf->GetHistNnn());
  fListBF->Add(bf->GetHistNpp());
  fListBF->Add(bf->GetHistNnp());
 
  //Balance function list: shuffling
  TList *fListBFS = 0x0;
  if(gRunShuffling) {
    fListBFS = new TList();
    fListBFS->SetName("listBFShuffled");
    fListBFS->SetOwner();
    fListBFS->Add(bfs->GetHistNp());
    fListBFS->Add(bfs->GetHistNn());
    fListBFS->Add(bfs->GetHistNpn());
    fListBFS->Add(bfs->GetHistNnn());
    fListBFS->Add(bfs->GetHistNpp());
    fListBFS->Add(bfs->GetHistNnp());
  }

  //Balance function list: event mixing
  TList *fListBFM = 0x0;
  if(gRunMixing) {
    fListBFM = new TList();
    fListBFM->SetName("listBFMixed");
    fListBFM->SetOwner();
    fListBFM->Add(bfm->GetHistNp());
    fListBFM->Add(bfm->GetHistNn());
    fListBFM->Add(bfm->GetHistNpn());
    fListBFM->Add(bfm->GetHistNnn());
    fListBFM->Add(bfm->GetHistNpp());
    fListBFM->Add(bfm->GetHistNnp());
  }
  //========================================================//

  //========================================================//
  //Event Mixing
  if(gRunMixing){
    Int_t trackDepth = 50000;
    Int_t poolsize   = 1000;  // Maximum number of events, ignored in the present implemented of AliEventPoolManager
        
    // centrality bins
    Double_t centralityBins[] = {0.,100.};
    Double_t* centbins        = centralityBins;
    Int_t nCentralityBins     = sizeof(centralityBins) / sizeof(Double_t) - 1;
    
    // Zvtx bins
    Double_t vertexBins[] = {-10., 10.}; 
    Double_t* vtxbins     = vertexBins;
    Int_t nVertexBins     = sizeof(vertexBins) / sizeof(Double_t) - 1;
    
    AliEventPoolManager *fPoolMgr = 0x0;
    fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins);
  }
  //========================================================//

  //========================================================//
  //Create the TChain object: events
  TChain *eventChain = new TChain("events");
  StructEvent tStructEvents;
  eventChain->SetBranchAddress("event",&tStructEvents);
  //========================================================//

  //========================================================//
  //Create the TChain object: particles 
  TChain *particleChain = new TChain("particles");
  StructParticle tStructParticles;
  particleChain->SetBranchAddress("particle",&tStructParticles);
  //========================================================//

  //========================================================//
  //Fill the TChain with the files in the directory
  for(Int_t iFile = 1; iFile <= aEventFiles; iFile++) {
    TString filename = aEventDir.Data();
    filename += "/event00"; filename += iFile;
    filename += ".root";
    cout<<"Adding file "<<filename.Data()<<" to the chain..."<<endl;
    eventChain->Add(filename.Data());
    particleChain->Add(filename.Data());
  }
  //========================================================//
 
  //========================================================//
  //QA histograms
  //Event stats.
  TString gCutName[5] = {"Total","Offline trigger",
                         "Vertex","Analyzed","sel. Centrality"};
  TH1F *fHistEventStats = new TH1F("fHistEventStats",
				   "Event statistics;;N_{events}",
				   5,0.5,5.5);
  for(Int_t i = 1; i <= 5; i++)
    fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());
  fList->Add(fHistEventStats);

  //Number of accepted particles
  TH1F *fHistNumberOfAcceptedTracks = new TH1F("fHistNumberOfAcceptedTracks",";N_{acc.};Entries",4001,-0.5,4000.5);
  fList->Add(fHistNumberOfAcceptedTracks);

  //Particle level: eta-phi
  TH2F *fHistEtaPhiPositive  = new TH2F("fHistEtaPhiPositive",";#eta;#varphi (rad)",80,gEtaMin,gEtaMax,72,0.0,2.*TMath::Pi());
  fList->Add(fHistEtaPhiPositive);
  TH2F *fHistEtaPhiNegative  = new TH2F("fHistEtaPhiNegative",";#eta;#varphi (rad)",80,gEtaMin,gEtaMax,72,0.0,2.*TMath::Pi());
  fList->Add(fHistEtaPhiNegative);

  //Particle level: pt
  TH1F *fHistPtPositive  = new TH1F("fHistPtPositive",";p_{T} (GeV/c)",100,0.01,30.01);
  fList->Add(fHistPtPositive);
  TH1F *fHistPtNegative  = new TH1F("fHistPtNegative",";p_{T} (GeV/c)",100,0.01,30.01);
  fList->Add(fHistPtNegative);
  //========================================================//

  //========================================================//
  //loop over the events
  Int_t nEvents = eventChain->GetEntries();
  cout<<"========================================="<<endl;
  cout<<"Number of events in the chain: "<<nEvents<<endl;
  cout<<"========================================="<<endl;
  Int_t iParticleCounter = 0;
  Int_t nTotalParticles = 0;
  
  for(Int_t iEvent = 0; iEvent < nEvents; iEvent++) {
    //for(Int_t iEvent = 0; iEvent < 1; iEvent++) {
    eventChain->GetEntry(iEvent);

    //========================================================//
    //Create the TObjArray object to store the particles
    TObjArray* tracksAccepted = new TObjArray;
    tracksAccepted->SetOwner(kTRUE);
    //========================================================//

    Int_t nParticles = tStructEvents.entries;
    if(((iEvent+1)%100)==0)
      cout<<"Event: "<<iEvent+1<<" - ID: "<<tStructEvents.eventID<<" - Entries: "<<tStructEvents.entries<<endl;

    //========================================================//
    //Filling event level histos
    fHistEventStats->Fill(1);
    fHistEventStats->Fill(2);
    fHistEventStats->Fill(3);
    fHistEventStats->Fill(4);
    fHistEventStats->Fill(5);

    Int_t gNumberOfAcceptedParticles = 0;
    Double_t gReactionPlane = 0.;
    Double_t gCharge = 0.;
    //========================================================//
    //loop over particles
    for(Int_t iParticle = 0; iParticle < nParticles; iParticle++) {
      particleChain->GetEntry(nTotalParticles+iParticle);
      
      //========================================================//
      //consider only primordial particles
      if(!IsPhysicalPrimary(tStructParticles.pid,
			    tStructParticles.fathereid,
			    gCharge)) continue;

      //========================================================//
      //Calculate kinematic variables
      Double_t gPt = TMath::Sqrt(TMath::Power(tStructParticles.px,2) + 
				 TMath::Power(tStructParticles.py,2));
      Double_t gP = TMath::Sqrt(TMath::Power(tStructParticles.px,2) + 
				TMath::Power(tStructParticles.py,2) + 
				TMath::Power(tStructParticles.pz,2) );
      Double_t gEta = -100.;
      if(gP != tStructParticles.pz)
	gEta = 0.5*TMath::Log((gP + tStructParticles.pz)/(gP - tStructParticles.pz));
      Double_t gPhi = TMath::Pi()+TMath::ATan2(-tStructParticles.py,-tStructParticles.px);

      //========================================================//
      //Fill QA
      if(gCharge > 0) {
	fHistEtaPhiPositive->Fill(gEta,gPhi);
	fHistPtPositive->Fill(gPt);
      }
      else if(gCharge < 0) {
	fHistEtaPhiNegative->Fill(gEta,gPhi);
	fHistPtNegative->Fill(gPt);
      }
      //========================================================//

      //========================================================//
      //Apply cuts
      if((gEta > gEtaMax)||(gEta < gEtaMin)) continue;
      if((gPt > gPtMax)||(gPt < gPtMin)) continue;
      
      tracksAccepted->Add(new AliBFBasicParticle(gEta,gPhi,gPt,gCharge, 1.));
      gNumberOfAcceptedParticles += 1;
      
      iParticleCounter += 1;
      //cout<<"\t Particle counter: "<<iParticleCounter<<" - Particle in the event: "<<iParticle+1<<" - eventID: "<<tStructParticles.eventid<<" - pid: "<<tStructParticles.pid<<" - fatherpid: "<<tStructParticles.fatherpid<<" - rootpid: "<<tStructParticles.rootpid<<" - fathereid: "<<tStructParticles.fathereid<<" - eid: "<<tStructParticles.eid<<endl;
    }//particle loop
    
    //========================================================//
    // Event mixing (borrowed code from the task) 
    if (gRunMixing) {
      Int_t fMixingTracks = 50000;
      AliEventPool* pool = fPoolMgr->GetEventPool(1.,0.,0.);
      
      if (!pool) {
	AliFatal(Form("No pool found"));
      }
      else {
	//pool->SetDebug(1);
	if (pool->IsReady() || pool->NTracksInPool() > fMixingTracks / 10 || pool->GetCurrentNEvents() >= 5){ 
	  
	  Int_t nMix = pool->GetCurrentNEvents();
	  //cout << "nMix = " << nMix << " tracks in pool = " << pool->NTracksInPool() << endl;
	  	  
	  // Fill mixed-event histos here  
	  for (Int_t jMix=0; jMix<nMix; jMix++) {
	    TObjArray* tracksMixed = pool->GetEvent(jMix);
	    bfm->CalculateBalance(gReactionPlane,tracksAccepted,tracksMixed,1.,0.);
	  }
	}
      
	// Update the Event pool
	pool->UpdatePool(tracksAccepted);
	//pool->PrintInfo();
	
      }//pool NULL check  
    }//run mixing
      //========================================================//

    //========================================================//
    // calculate balance function
    bf->CalculateBalance(gReactionPlane,tracksAccepted,NULL,1.,0.);

    fHistNumberOfAcceptedTracks->Fill(gNumberOfAcceptedParticles);
    nTotalParticles += nParticles;
  }

  //========================================================//
  //Output file
  TFile *f = TFile::Open("AnalysisResults.root","recreate");
   TDirectoryFile *dir = new TDirectoryFile("PWGCFEbyE.outputBalanceFunctionPsiAnalysis","PWGCFEbyE.outputBalanceFunctionPsiAnalysis");
  
  fList->SetName("listQA");
  fList->SetOwner(kTRUE);
  dir->Add(fList);

  fListBF->SetName("listBF");
  fListBF->SetOwner(kTRUE);
  dir->Add(fListBF);

  if(gRunMixing) {
    fListBFM->SetName("listBFMixed");
    fListBFM->SetOwner(kTRUE);
    dir->Add(fListBFM);
  }
  
  dir->Write(dir->GetName(),TObject::kSingleKey);
  f->Close();  
  //========================================================//

  // Print real and CPU time used for analysis:  
  timer.Stop();
  timer.Print();
}

//______________________________________________________________//
Bool_t IsPhysicalPrimary(Int_t gPDGCode,
			 Int_t gFathereid,
			 Double_t &gCharge) {
  //Check whether the primordial particle belongs to the list 
  //of known particles
  Bool_t kStatus = kFALSE;
  if(gFathereid != -1) 
    return kStatus;

  //List of stable particles
  const Int_t kNstable = 3;
  Int_t pdgStable[kNstable] = {
    211,         // Pion
    321,         // Kaon
    2212,        // Proton 
  };
  
  if(gPDGCode < 0) gCharge = -1;
  else if(gPDGCode > 0) gCharge = 1;

  for(Int_t iParticle = 0; iParticle < kNstable; iParticle++) {
    if((TMath::Abs(gPDGCode) == pdgStable[iParticle])) {
      kStatus = kTRUE;
      return kStatus;
    }
  }

  return kFALSE;
}
Commit	Line	Data
ee97f33d	1	//========================================================//
	2	//Macro that reads the output of the hydro calculations
	3	//(Therminator) using the format that can be found at
	4	//arXiv:1102.0273
	5	//Author: Panos.Christakoglou@nikhef.nl
	6	//========================================================//
	7
	8	//========================================================//
	9	//Event structure
	10	struct StructEvent {
	11	UInt_t eventID;//unique identifier of the event (random number)
	12	UInt_t entries;//number of particles of each event
	13	UInt_t entriesprev;//set to 0 by default
	14	};
	15	//========================================================//
	16
	17	//========================================================//
	18	//Particle structure
a7f9d873	19	//Primordial particles: fathereid==-1 ==>
a7f9d873	20	//pid==fatherpid==rootpid
ee97f33d	21	struct StructParticle {
	22	Float_t mass;//the mass of the particle (in GeV)
	23	Float_t t;//the time coordinate of the particle in fm/c
	24	Float_t x;//the spacial coordinate x in fm/c
	25	Float_t y;//the spacial coordinate y in fm/c
	26	Float_t z;//the spacial coordinate z in fm/c
	27	Float_t e;//the energy in GeV
	28	Float_t px;//the x-coordinate of the particle's momentum in GeV
	29	Float_t py;//the y-coordinate of the particle's momentum in GeV
	30	Float_t pz;//the z-coordinate of the particle's momentum in GeV
	31	Int_t decayed;//decay flag (1: decayed, 0: no decay)
	32	Int_t pid;//PDG of particle
	33	Int_t fatherpid;//PDG of parent
	34	Int_t rootpid;//root (primordial) particle PDG number
	35	Int_t eid;//sequence number in the event
	36	Int_t fathereid;//parent sequence number in the event
	37	UInt_t eventid;//unique identifier of the event (random number)
	38	};
	39	//========================================================//
	40
	41	//========================================================//
	42	//Balance function analysis variables
	43	Bool_t gRunShuffling=kFALSE;
eaae9cd6	44	Bool_t gRunMixing=kTRUE;
ee97f33d	45	Bool_t gRunMixingWithEventPlane=kFALSE;
	46
	47	Double_t gEtaMin = -0.8;
	48	Double_t gEtaMax = 0.8;
	49	Double_t gPtMin = 0.2;
	50	Double_t gPtMax = 20.0;
	51	//========================================================//
	52
	53
8e81f4e0	54	void runBalanceFunctionOnHydro(TString aEventDir = ".",
8e81f4e0	55	Int_t aEventFiles = 9) {
ee97f33d	56	//Macro that reads the themrinator events
	57	//Author: Panos.Christakoglou@nikhef.nl
	58	// Time:
	59	TStopwatch timer;
	60	timer.Start();
	61
	62	//========================================================//
	63	//Load the aliroot libraries
4070f709	64	gSystem->Load("libSTEERBase");
	65	gSystem->Load("libESD");
	66	gSystem->Load("libAOD");
	67	gSystem->Load("libANALYSIS");
	68	gSystem->Load("libANALYSISalice");
	69	gSystem->Load("libEventMixing");
	70	gSystem->Load("libCORRFW");
	71	gSystem->Load("libPWGTools");
	72	gSystem->Load("libPWGCFebye");
ee97f33d	73	//========================================================//
	74
	75	//========================================================//
	76	//Configure the bf objects
	77	AliBalancePsi *bf = new AliBalancePsi();
	78	bf->SetAnalysisLevel("MC");
	79	bf->SetShuffle(gRunShuffling);
	80	bf->SetEventClass("EventPlane");
	81	bf->SetDeltaEtaMax(TMath::Abs(gEtaMax-gEtaMin));
	82	bf->InitHistograms();
	83
	84	AliBalancePsi *bfs = 0x0;
	85	if(gRunShuffling) {
	86	bfs = new AliBalancePsi();
	87	bfs->SetAnalysisLevel("MC");
	88	bfs->SetEventClass("EventPlane");
	89	bfs->SetDeltaEtaMax(TMath::Abs(gEtaMax-gEtaMin));
	90	bfs->InitHistograms();
	91	}
	92
	93	AliBalancePsi *bfm = 0x0;
	94	if(gRunMixing) {
	95	bfm = new AliBalancePsi();
	96	bfm->SetAnalysisLevel("MC");
	97	bfm->SetShuffle(gRunShuffling);
	98	bfm->SetEventClass("EventPlane");
	99	bfm->SetDeltaEtaMax(TMath::Abs(gEtaMax-gEtaMin));
	100	bfm->InitHistograms();
	101	}
	102	//========================================================//
	103
	104	//========================================================//
	105	//Output objects
	106	//QA list
	107	fList = new TList();
	108	fList->SetName("listQA");
	109	fList->SetOwner();
	110
	111	//Balance Function list
	112	TList *fListBF = new TList();
	113	fListBF->SetName("listBF");
	114	fListBF->SetOwner();
eaae9cd6	115	fListBF->Add(bf->GetHistNp());
	116	fListBF->Add(bf->GetHistNn());
	117	fListBF->Add(bf->GetHistNpn());
	118	fListBF->Add(bf->GetHistNnn());
	119	fListBF->Add(bf->GetHistNpp());
	120	fListBF->Add(bf->GetHistNnp());
	121
ee97f33d	122	//Balance function list: shuffling
	123	TList *fListBFS = 0x0;
	124	if(gRunShuffling) {
	125	fListBFS = new TList();
	126	fListBFS->SetName("listBFShuffled");
	127	fListBFS->SetOwner();
	128	fListBFS->Add(bfs->GetHistNp());
	129	fListBFS->Add(bfs->GetHistNn());
	130	fListBFS->Add(bfs->GetHistNpn());
	131	fListBFS->Add(bfs->GetHistNnn());
	132	fListBFS->Add(bfs->GetHistNpp());
	133	fListBFS->Add(bfs->GetHistNnp());
	134	}
	135
	136	//Balance function list: event mixing
	137	TList *fListBFM = 0x0;
	138	if(gRunMixing) {
	139	fListBFM = new TList();
	140	fListBFM->SetName("listBFMixed");
	141	fListBFM->SetOwner();
	142	fListBFM->Add(bfm->GetHistNp());
	143	fListBFM->Add(bfm->GetHistNn());
	144	fListBFM->Add(bfm->GetHistNpn());
	145	fListBFM->Add(bfm->GetHistNnn());
	146	fListBFM->Add(bfm->GetHistNpp());
	147	fListBFM->Add(bfm->GetHistNnp());
	148	}
	149	//========================================================//
	150
	151	//========================================================//
	152	//Event Mixing
	153	if(gRunMixing){
	154	Int_t trackDepth = 50000;
	155	Int_t poolsize = 1000; // Maximum number of events, ignored in the present implemented of AliEventPoolManager
eaae9cd6	156
	157	// centrality bins
	158	Double_t centralityBins[] = {0.,100.};
	159	Double_t* centbins = centralityBins;
	160	Int_t nCentralityBins = sizeof(centralityBins) / sizeof(Double_t) - 1;
	161
	162	// Zvtx bins
	163	Double_t vertexBins[] = {-10., 10.};
	164	Double_t* vtxbins = vertexBins;
	165	Int_t nVertexBins = sizeof(vertexBins) / sizeof(Double_t) - 1;
ee97f33d	166
ee97f33d	167	AliEventPoolManager *fPoolMgr = 0x0;
eaae9cd6	168	fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins);
ee97f33d	169	}
	170	//========================================================//
	171
	172	//========================================================//
	173	//Create the TChain object: events
	174	TChain *eventChain = new TChain("events");
	175	StructEvent tStructEvents;
	176	eventChain->SetBranchAddress("event",&tStructEvents);
	177	//========================================================//
	178
	179	//========================================================//
	180	//Create the TChain object: particles
	181	TChain *particleChain = new TChain("particles");
	182	StructParticle tStructParticles;
	183	particleChain->SetBranchAddress("particle",&tStructParticles);
	184	//========================================================//
	185
	186	//========================================================//
	187	//Fill the TChain with the files in the directory
8e81f4e0	188	for(Int_t iFile = 1; iFile <= aEventFiles; iFile++) {
ee97f33d	189	TString filename = aEventDir.Data();
	190	filename += "/event00"; filename += iFile;
	191	filename += ".root";
	192	cout<<"Adding file "<<filename.Data()<<" to the chain..."<<endl;
	193	eventChain->Add(filename.Data());
	194	particleChain->Add(filename.Data());
	195	}
	196	//========================================================//
	197
	198	//========================================================//
eaae9cd6	199	//QA histograms
ee97f33d	200	//Event stats.
	201	TString gCutName[5] = {"Total","Offline trigger",
	202	"Vertex","Analyzed","sel. Centrality"};
a7f9d873	203	TH1F *fHistEventStats = new TH1F("fHistEventStats",
	204	"Event statistics;;N_{events}",
	205	5,0.5,5.5);
ee97f33d	206	for(Int_t i = 1; i <= 5; i++)
	207	fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());
	208	fList->Add(fHistEventStats);
	209
	210	//Number of accepted particles
a7f9d873	211	TH1F *fHistNumberOfAcceptedTracks = new TH1F("fHistNumberOfAcceptedTracks",";N_{acc.};Entries",4001,-0.5,4000.5);
ee97f33d	212	fList->Add(fHistNumberOfAcceptedTracks);
eaae9cd6	213
	214	//Particle level: eta-phi
	215	TH2F fHistEtaPhiPositive = new TH2F("fHistEtaPhiPositive",";#eta;#varphi (rad)",80,gEtaMin,gEtaMax,72,0.0,2.TMath::Pi());
	216	fList->Add(fHistEtaPhiPositive);
	217	TH2F fHistEtaPhiNegative = new TH2F("fHistEtaPhiNegative",";#eta;#varphi (rad)",80,gEtaMin,gEtaMax,72,0.0,2.TMath::Pi());
	218	fList->Add(fHistEtaPhiNegative);
	219
	220	//Particle level: pt
	221	TH1F *fHistPtPositive = new TH1F("fHistPtPositive",";p_{T} (GeV/c)",100,0.01,30.01);
	222	fList->Add(fHistPtPositive);
	223	TH1F *fHistPtNegative = new TH1F("fHistPtNegative",";p_{T} (GeV/c)",100,0.01,30.01);
	224	fList->Add(fHistPtNegative);
ee97f33d	225	//========================================================//
	226
	227	//========================================================//
	228	//loop over the events
	229	Int_t nEvents = eventChain->GetEntries();
	230	cout<<"========================================="<<endl;
	231	cout<<"Number of events in the chain: "<<nEvents<<endl;
	232	cout<<"========================================="<<endl;
	233	Int_t iParticleCounter = 0;
	234	Int_t nTotalParticles = 0;
	235
eaae9cd6	236	for(Int_t iEvent = 0; iEvent < nEvents; iEvent++) {
eaae9cd6	237	//for(Int_t iEvent = 0; iEvent < 1; iEvent++) {
ee97f33d	238	eventChain->GetEntry(iEvent);
ee97f33d	239
a7f9d873	240	//========================================================//
	241	//Create the TObjArray object to store the particles
	242	TObjArray* tracksAccepted = new TObjArray;
	243	tracksAccepted->SetOwner(kTRUE);
	244	//========================================================//
	245
ee97f33d	246	Int_t nParticles = tStructEvents.entries;
a7f9d873	247	if(((iEvent+1)%100)==0)
a7f9d873	248	cout<<"Event: "<<iEvent+1<<" - ID: "<<tStructEvents.eventID<<" - Entries: "<<tStructEvents.entries<<endl;
ee97f33d	249
a7f9d873	250	//========================================================//
	251	//Filling event level histos
	252	fHistEventStats->Fill(1);
	253	fHistEventStats->Fill(2);
	254	fHistEventStats->Fill(3);
	255	fHistEventStats->Fill(4);
	256	fHistEventStats->Fill(5);
	257
	258	Int_t gNumberOfAcceptedParticles = 0;
	259	Double_t gReactionPlane = 0.;
eaae9cd6	260	Double_t gCharge = 0.;
ee97f33d	261	//========================================================//
	262	//loop over particles
	263	for(Int_t iParticle = 0; iParticle < nParticles; iParticle++) {
	264	particleChain->GetEntry(nTotalParticles+iParticle);
	265
a7f9d873	266	//========================================================//
a7f9d873	267	//consider only primordial particles
eaae9cd6	268	if(!IsPhysicalPrimary(tStructParticles.pid,
	269	tStructParticles.fathereid,
	270	gCharge)) continue;
a7f9d873	271
	272	//========================================================//
	273	//Calculate kinematic variables
	274	Double_t gPt = TMath::Sqrt(TMath::Power(tStructParticles.px,2) +
	275	TMath::Power(tStructParticles.py,2));
	276	Double_t gP = TMath::Sqrt(TMath::Power(tStructParticles.px,2) +
	277	TMath::Power(tStructParticles.py,2) +
	278	TMath::Power(tStructParticles.pz,2) );
	279	Double_t gEta = -100.;
	280	if(gP != tStructParticles.pz)
	281	gEta = 0.5*TMath::Log((gP + tStructParticles.pz)/(gP - tStructParticles.pz));
eaae9cd6	282	Double_t gPhi = TMath::Pi()+TMath::ATan2(-tStructParticles.py,-tStructParticles.px);
	283
	284	//========================================================//
	285	//Fill QA
	286	if(gCharge > 0) {
	287	fHistEtaPhiPositive->Fill(gEta,gPhi);
	288	fHistPtPositive->Fill(gPt);
	289	}
	290	else if(gCharge < 0) {
	291	fHistEtaPhiNegative->Fill(gEta,gPhi);
	292	fHistPtNegative->Fill(gPt);
	293	}
	294	//========================================================//
	295
a7f9d873	296	//========================================================//
	297	//Apply cuts
	298	if((gEta > gEtaMax)\|\|(gEta < gEtaMin)) continue;
	299	if((gPt > gPtMax)\|\|(gPt < gPtMin)) continue;
	300
	301	tracksAccepted->Add(new AliBFBasicParticle(gEta,gPhi,gPt,gCharge, 1.));
	302	gNumberOfAcceptedParticles += 1;
	303
ee97f33d	304	iParticleCounter += 1;
a7f9d873	305	//cout<<"\t Particle counter: "<<iParticleCounter<<" - Particle in the event: "<<iParticle+1<<" - eventID: "<<tStructParticles.eventid<<" - pid: "<<tStructParticles.pid<<" - fatherpid: "<<tStructParticles.fatherpid<<" - rootpid: "<<tStructParticles.rootpid<<" - fathereid: "<<tStructParticles.fathereid<<" - eid: "<<tStructParticles.eid<<endl;
ee97f33d	306	}//particle loop
a7f9d873	307
	308	//========================================================//
	309	// Event mixing (borrowed code from the task)
eaae9cd6	310	if (gRunMixing) {
	311	Int_t fMixingTracks = 50000;
	312	AliEventPool* pool = fPoolMgr->GetEventPool(1.,0.,0.);
a7f9d873	313
a7f9d873	314	if (!pool) {
eaae9cd6	315	AliFatal(Form("No pool found"));
a7f9d873	316	}
a7f9d873	317	else {
eaae9cd6	318	//pool->SetDebug(1);
	319	if (pool->IsReady() \|\| pool->NTracksInPool() > fMixingTracks / 10 \|\| pool->GetCurrentNEvents() >= 5){
	320
	321	Int_t nMix = pool->GetCurrentNEvents();
	322	//cout << "nMix = " << nMix << " tracks in pool = " << pool->NTracksInPool() << endl;
	323
	324	// Fill mixed-event histos here
	325	for (Int_t jMix=0; jMix<nMix; jMix++) {
	326	TObjArray* tracksMixed = pool->GetEvent(jMix);
	327	bfm->CalculateBalance(gReactionPlane,tracksAccepted,tracksMixed,1.,0.);
	328	}
	329	}
a7f9d873	330
eaae9cd6	331	// Update the Event pool
	332	pool->UpdatePool(tracksAccepted);
	333	//pool->PrintInfo();
	334
a7f9d873	335	}//pool NULL check
eaae9cd6	336	}//run mixing
a7f9d873	337	//========================================================//
	338
	339	//========================================================//
	340	// calculate balance function
eaae9cd6	341	bf->CalculateBalance(gReactionPlane,tracksAccepted,NULL,1.,0.);
a7f9d873	342
a7f9d873	343	fHistNumberOfAcceptedTracks->Fill(gNumberOfAcceptedParticles);
ee97f33d	344	nTotalParticles += nParticles;
	345	}
	346
	347	//========================================================//
	348	//Output file
a7f9d873	349	TFile *f = TFile::Open("AnalysisResults.root","recreate");
8e81f4e0	350	TDirectoryFile *dir = new TDirectoryFile("PWGCFEbyE.outputBalanceFunctionPsiAnalysis","PWGCFEbyE.outputBalanceFunctionPsiAnalysis");
	351
	352	fList->SetName("listQA");
	353	fList->SetOwner(kTRUE);
	354	dir->Add(fList);
	355
	356	fListBF->SetName("listBF");
	357	fListBF->SetOwner(kTRUE);
	358	dir->Add(fListBF);
	359
	360	if(gRunMixing) {
	361	fListBFM->SetName("listBFMixed");
	362	fListBFM->SetOwner(kTRUE);
	363	dir->Add(fListBFM);
	364	}
	365
	366	dir->Write(dir->GetName(),TObject::kSingleKey);
ee97f33d	367	f->Close();
	368	//========================================================//
	369
	370	// Print real and CPU time used for analysis:
	371	timer.Stop();
	372	timer.Print();
	373	}
eaae9cd6	374
	375	//______________________________________________________________//
	376	Bool_t IsPhysicalPrimary(Int_t gPDGCode,
	377	Int_t gFathereid,
	378	Double_t &gCharge) {
	379	//Check whether the primordial particle belongs to the list
	380	//of known particles
	381	Bool_t kStatus = kFALSE;
	382	if(gFathereid != -1)
	383	return kStatus;
	384
	385	//List of stable particles
	386	const Int_t kNstable = 3;
	387	Int_t pdgStable[kNstable] = {
	388	211, // Pion
	389	321, // Kaon
	390	2212, // Proton
	391	};
	392
	393	if(gPDGCode < 0) gCharge = -1;
	394	else if(gPDGCode > 0) gCharge = 1;
	395
	396	for(Int_t iParticle = 0; iParticle < kNstable; iParticle++) {
	397	if((TMath::Abs(gPDGCode) == pdgStable[iParticle])) {
	398	kStatus = kTRUE;
	399	return kStatus;
	400	}
	401	}
	402
	403	return kFALSE;
	404	}